Rosalie Plantefève
email: rosalie.plantefeve@inria.fr

Igor Peterlik

Nazim Haouchine

Stéphane Cotin

Patient-specific Biomechanical Modeling for Guidance during Minimally-invasive Hepatic Surgery

Keywords: patient-specific modeling, non-rigid registration, minimally invasive surgery, real-time simulation

During the minimally-invasive liver surgery, only the partial surface view of the liver is usually provided to the surgeon via the laparoscopic camera. Therefore, it is necessary to estimate the actual position of the internal structures such as tumors and vessels from the pre-operative images. Nevertheless, such task can be highly challenging since during the intervention, the abdominal organs undergo important deformations due to the pneumoperitoneum, respiratory and cardiac motion and the interaction with the surgical tools. Therefore, a reliable automatic system for intra-operative guidance requires fast and reliable registration of the pre-and intra-operative data.

In this paper we present a complete pipeline for the registration of pre-operative patient-specific image data to the sparse and incomplete intra-operative data. While the intra-operative data is represented by a point cloud extracted from the stereo-endoscopic images, the pre-operative data is used to reconstruct a biomechanical model which is necessary for accurate estimation of the position of the internal structures, considering the actual deformations. This model takes into account the patient-specific liver anatomy composed of parenchyma, vascularization and capsule, and is enriched with anatomical boundary conditions transferred from an atlas. The registration process employs the iterative closest point technique together with a penalty-based method.

We perform a quantitative assessment based on the evaluation of the target registration error on synthetic data as well as a qualitative assessment on real patient data. We demonstrate that the proposed registration method provides good results in terms of both accuracy and robustness w. r. t. the quality of the intra-operative data.

Introduction

In the last decades, advances in medicine have seen the emergence of minimally invasive surgery. In this surgical approach, the surgeon manipulates the organs with instruments inserted through trocars placed in small abdominal incisions. The view of the operating field is provided by a laparoscopic camera, inserted also through a trocar, allowing the surgeon to perform the surgery by watching a high-definition monitor placed above the patient. Minimally invasive surgery provides real benefits to the patient, by reducing pain, bleeding and risks of infection and therefore shortening recovery time. However, since the manipulation of instruments is indirect and the visual feedback is limited, it remains quite challenging from a surgical standpoint. Improved navigation, planning, and visualization of internal structures are among the main requirements. The use of a laparoscopic camera in the clinical routine has naturally led the research community to investigate registration techniques in order to overlay the pre-operative anatomical models directly onto the intra-operative view. In the context of hepatic surgery, the objective is to locate tumors with greater accuracy and preserve vessels which are needed for the post-operative regeneration of the liver tissue. For this reason, and also due to the important deformations underwent by the liver during the procedure, minimally invasive hepatic surgery is a challenge for intra-operative guidance methods.

Registration between intra-operative data and pre-operative data is a cross-disciplinary research domain combining disciplines such as computer vision, image processing, computer graphics and computational mechanics methods. Indeed, where surgical vision methods [START_REF] Stoyanov | Surgical vision[END_REF] help process the endoscopic images in order to recover organs shapes, perform instruments tracking, or estimate camera position, computer graphics techniques aim at providing three-dimensional virtual representations of real organs using advanced geometri-Fig. 1: Overall concept of our approach illustrating the non-rigid registration process using a patient-specific biomechanical model derived from pre-operative data. The registration with the endoscopic image permits the visualization of internal structures during the operation in order to provide surgeons with decision support. cal modeling and rendering algorithms. These three-dimensional representations are often considered as purely geometric, since most existing surgical navigation systems assume rigid organ motion [START_REF] Nicolau | Augmented reality in laparoscopic surgical oncology[END_REF]. Nevertheless, in real clinical routines, this assumption no longer holds, given the elastic characteristics of soft tissues. In addition, in minimally invasive surgery, the abdomen is inflated with CO 2 , creating large deformations of the organs to which respiratory and cardiac motions are added, as well as interactions with surgical instruments. As a consequence, more advanced modeling methods need to be envisioned, such that elastic deformations linking pre-and intra-operative data can be simulated in real-time.

The benefits of surgery navigation are particularly clear in laparoscopic hepatectomy where it can enable the visualization of pre-operative data directly fused with endoscopic images. Generally, a geometrical model constructed from pre-operative data is used during the navigation [START_REF] Puerto-Souza | Toward long-term and accurate augmented-reality display for minimally-invasive surgery[END_REF][START_REF] Su | Augmented reality during robot-assisted laparoscopic partial nephrectomy: Toward real-time 3d-ct to stereoscopic video registration[END_REF][START_REF] Nicolau | Augmented reality in laparoscopic surgical oncology[END_REF]. Although such approaches offer the possibility to achieve realtime registration and give visually coherent augmentation, they are limited to surface overlay, as neither the vessels nor the tumors are considered. They also rarely account for soft tissue deformation. To provide in-depth structures overlay, elastic registration and volume computation has to be considered [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF][START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF][START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF]. Physically-based elastic models have proven to be relevant for volume deformation tracking since they permit to estimate the in-depth motion from surface deformation. In [START_REF] Pratt | Dynamic guidance for robotic surgery using image-constrained biomechanical models[END_REF], a 4D scan of the heart is coupled with a biomechanical model. It is controlled by surface constraints created by features extracted from a stereo-vision camera stream and allows for a quite accurate estimation of the position of internal structures. This approach is, however, dedicated to cyclic movements and cannot be translated to liver surgery where large and unpredictable deformations may occur. Based on a dynamic model, Schaerer et al. [START_REF] Schaerer | A dynamic elastic model for segmentation and tracking of the heart in mr image sequences[END_REF] proposed to use forces measured in the image to drive the model towards the boundaries of the object. The main limitation of this method is that it assumes the organ to be homogeneous, and that the entire organ is visible in the MRI data. A Neohookean elastic model is employed in [START_REF] Speidel | Intraoperative surface reconstruction and biomechanical modeling for soft tissue registration[END_REF] in order to perform intra-operative registration between stereo-endoscopic images and pre-operative model of the liver. The accuracy of registration is assessed on a phantom and shows the suitability of the model for an accurate registration, nevertheless not suitable for a real time simulation. In our previous work [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF], we used a deformable biomechanical model accounting for heterogeneity and anisotropy for hepatic surgery guidance. Convincing results are reported on both in vivo and phantom data under well-defined conditions. Nevertheless, the boundary conditions, which are essential to yield good re-sults, were defined in an ad-hoc manner during the initialization phase. In [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF], shape-matching is used for the registration of a physically-based model derived from pre-operative mesh and intra-operative surface computed from camera. The intra-operative organ shape is modeled as an electrostatic elastic problem. The elastic model is charged electrically to slide into an oppositely charged organ shape representation. Nonetheless, the method uses a simplified homogeneous model of the liver and it is assumed that at least 50% of the organ surface is captured by the camera, which is generally not possible. A pose-independent matching of intra-operative data acquired by a camera was proposed in [START_REF] Santos | Pose-independent surface matching for intraoperative soft-tissue marker-less registration[END_REF]. The method is automatic and takes into account data noise and tissue deformations. However, the approach requires to have access to an advanced camera hardware providing depth information, which facilitates the reconstruction of the surface, but is not usable in clinical (minimally invasive) conditions. Based on the work of [START_REF] Bano | Simulation of pneumoperitoneum for laparoscopic surgery planning[END_REF], Oktay et al. [START_REF] Oktay | Biomechanically driven registration of pre-to intra-operative 3d images for laparoscopic surgery[END_REF] proposed to combine the intra-operative data acquired from CT-scans with a biomechanical simulation of gas insufflation to accurately perform the registration. Although this method provides accurate registration, it relies on intra-operative scans, which are currently not available during clinical routines. Clements et al. [START_REF] Clements | Robust surface registration using salient anatomical features for image-guided liver surgery: Algorithm and validation[END_REF] introduced an approach that exploits salient anatomical features, identifiable in both the pre-operative images and intra-operative liver surface. Based on an Iterative Closest Point (ICP) the method performs well, but is limited to rigid transformations. In [START_REF] Plantefeve | Automatic of pre and intraoperative data using anatomical landmarks for augmented laparoscopic liver surgery[END_REF] an automatic registration of the liver using only the endoscopic images is proposed. It relies on an approach similar to [START_REF] Clements | Robust surface registration using salient anatomical features for image-guided liver surgery: Algorithm and validation[END_REF] but using a deformable model. Nonetheless, the method performs the registration only with a static point cloud, thus not handling liver motion due to respiratory and cardiac motions.

While the above approaches are often based on a biomechanical model to describe the behavior of the deformable body, the modeling of the organ remains relatively simple and usually not compatible with the real-time simulation. For instance, in the case of the liver, it has been demonstrated experimentally that this organ has a complex anatomy with a significant level of patient-specificity [START_REF] Yeh | Elastic modulus measurements of human liver and correlation with pathology[END_REF][START_REF] Mazza | The mechanical response of human liver and its relation to histology: An in vivo study[END_REF]. In our context, using a more accurate model will improve the accuracy of the internal structure position estimation [START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF]. Such an improvement in the model can be obtained by accounting for the three different constituents (parenchyma, vessels, Glisson's capsule) of the liver, each having its characteristic mechanical properties [START_REF] Umale | Experimental in vitro mechanical characterization of porcine Glisson's capsule and hepatic veins[END_REF]. Another element influencing the behavior of the liver is the pathology, which in certain cases can stiffen the organ significantly. Among the important body of work that exists regarding the biomechanical properties of the liver parenchyma, only few studies focused on the role of vascularization inside the tissue. In [START_REF] Kerdok | Effects of perfusion on the viscoelastic characteristics of liver[END_REF], a visco-elastic model of the liver is proposed as well as mate-rial parameters experimentally measured ex vivo on perfused liver. In [START_REF] Nguyen | Patient Specific Biomechanical Modeling of Hepatic Vasculature for Augmented Reality Surgery[END_REF], a patient specific model of hepatic vasculature is proposed. The material properties of vessels are modelled by non-linear constitutive law. Nevertheless, the model does not allow for real-time performance as the vessel walls are modelled with large number of finite elements.

This article is an extension of our previous works [START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF][START_REF] Plantefeve | Automatic of pre and intraoperative data using anatomical landmarks for augmented laparoscopic liver surgery[END_REF], where we improve upon [START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF] by adding anatomically correct boundary conditions, and extend the work presented in [START_REF] Plantefeve | Automatic of pre and intraoperative data using anatomical landmarks for augmented laparoscopic liver surgery[END_REF] by considering a non-static point cloud in the registration process. In this paper we mainly focus on the usability and assessment of the framework depicted in Fig. 1. In particular we emphasize the patient-specific model of the liver, which takes into account the vascular network as the main source of heterogeneity and anisotropy, yet allows for real-time augmentation and is directly derived from medical images routinely acquired in the case of liver cancer. The details about the model reconstruction are presented in section 2.1. The biomechanical model is constrained by anatomically correct boundary conditions, which are transferred from a statistical atlas. These boundary conditions permit to obtain a well-constrained system during the registration process which is based on a two-step approach. First, an iterative algorithm brings the pre-operative model in the camera reference frame and deforms it according to the point cloud shape. Then, a tracking algorithm guarantees the real-time registration of the biomechanical model with the liver view of the endoscopic video flow. The method of both the initial and the real-time registrations are described in section 2.2. Further, the section 3 provides an assessment of the method using the synthetic and in vivo experimental data. The limitations of the method as well as its applicability are finally discussed in section 4.

Materials and Methods

Patient-specific Composite Model of Human Liver

The biomechanical model used for the initial registration and intra-operative positional prediction of the internal structures (such as vessels and lesions) has two aspects: the reconstruction of the geometry and the physical formulation of the problem. It has been demonstrated via rheological experiments presented in [START_REF] Umale | Experimental in vitro mechanical characterization of porcine Glisson's capsule and hepatic veins[END_REF] that the biomechanical behaviour of the liver is determined by three constituents: parenchyma, Glisson's capsule and vascularization. Therefore, our framework relies on the composite model presented in [START_REF] Peterlik | Modeling and real-time simulation of a vascularized liver tissue[END_REF][START_REF] Peterlík | Complete real-time liver model including glissons capsule, vascularization and parenchyma[END_REF] which combines the mechanical response of the three constituents while allows for real-time employment of the resulting virtual object.

In this section, the semi-automatic pipeline which creates the geometrical representation from the input images (such as pre-operative CT or MRI) is first described. Further, the mechanical aspects of the composite model are presented together with an efficient solution method.

Reconstruction of the model geometry

The geometry of the model is obtained by a semi-automatic pipeline presented in Fig. 2. Starting from the pre-operative images, the pipeline has two branches: the geometrical reconstruction of the parenchyma and the capsule and the reconstruction of the vascular model. While the former consists of segmentation and meshing, the latter requires skeletonization of the vascular tree and construction of the geometrical representation suitable for the finite-element formulation based on beam elements. Below we present details for each phase of the geometrical reconstruction.

Geometrical reconstruction of parenchyma and capsule

To perform the liver segmentation, we employ the field level set algorithm presented in [START_REF] Yushkevich | User-guided 3d active contour segmentation of anatomical structures: significantly improved efficiency and reliability[END_REF]. The quality of the segmented map is improved by several filtering methods, namely by Laplacian smoothing, hole-filling and detection of the continuous components. Afterwards, it is used to generate the finite element mesh for both the parenchyma and the capsule. In our pipeline, we employ the method presented in [START_REF] Boltcheva | Mesh generation from 3d multi-material images[END_REF] which generates the 3D tetrahedral mesh of the parenchyma directly from the segmented map and reduces the number of degenerated elements and slivers using mesh optimizations such as Lloyd and ODT smoothers. The surface triangulation of the mesh is extracted from the tetrahedral mesh, resulting in the geometrical domain of the capsule.

Geometrical reconstruction of the vascularization

Blood vessels are segmented from the same image modality as the parenchyma by a semi-automatic method [START_REF] Yushkevich | User-guided 3d active contour segmentation of anatomical structures: significantly improved efficiency and reliability[END_REF]. Then a skeletonization is computed using an algorithm based on Voronoi diagram [START_REF] Antiga | Centerline computation and geometric analysis of branching tubular surfaces with application to blood vessel modeling[END_REF] which as the input takes the surface of the segmented map. This surface is usually constructed using a marching cubes algorithm, however the resulting skeleton often suffers from incompleteness as the method is very sensitive to the quality of the segmentation. Therefore, we adapted the method based on iterative Dijkstra minimum cost spanning tree presented in [START_REF] Verscheure | Three-dimensional skeletonization and symbolic description in vascular imaging: preliminary results[END_REF]. The method works directly with the segmented map which is converted into a weighted 3D graph where each voxel belonging to the segmented volume is prepresented by a weighted node.

The method constructs the skeleton in two steps: first, the spanning tree is constructed iteratively starting from the root voxel. The edges between voxels are constructed recursively using a sorted heap: in each step, the head of the heap having the minimal weight is marked and all its unmarked neighbors are inserted into the heap. The node weight used to sort the heap is defined as 1 r b where r b is a distance from boundary which is obtained from distance map pre-computed from the binary image of the vascular tree. The weight thus prioritizes the voxels with higher probability to lie on the centerline, i. e. further from the vessel boundary. Moreover, for each visited voxel, its distance from the root (d r ) is stored, being the geometrical distance along the edges contructed so far from the root.

Second, the centerlines are extracted recursively from the tree:

-Path P 0 of order 0 is constructed as the shortest path connecting the root and the voxel with the highest value of d r in the graph. -Path P n of order n is constructed in two steps:

voxel expansion: for each voxel v of the path P n-1 of order (n-1), a set V v of all voxels accessible from v is constructed. path extraction: the shortest path (v, w) is constructed:

v is the parent voxel (on path of order (n -1)) and w is the voxel from V v having the highest value of d r .

We introduced several modifications w. r. t. the original method [START_REF] Verscheure | Three-dimensional skeletonization and symbolic description in vascular imaging: preliminary results[END_REF]: first of all, in order to increase the smoothness of the resulting centerline, we evaluated the extraction for different exponents e ∈ {1, 2, 3, 4, 5, 6} in the definition of the node weight 1

r e b
. We obtained the best results in terms of centerline smoothness and continuity using e = 4. Further, we added two heuristics in order to improve the quality of the skeleton. The first heuristic attempts to shorten the new path P n (v, w) by finding a parent voxel u which is closer to the terminal voxel w. Mathematically, we look for a voxel u ∈ P n-1 such that: u is adjacent to a voxel from P n and for the lengths of the paths, it holds:

|(u, w)| < |(v, w)|. If such voxel exists, path (v, w) is replaced by path (u, w).
The second heuristic detects the false paths which lie inside a tube in which a path already exists. Before testing the new path P n (v, w), we find a voxel u ∈ P n-1 which lies in the parent path and minimizes the Euclidean distance |(u, w)|. This distance is not measured along the path but it is simply the length of a straight line connecting u and w. A new path P n (v, w) is removed if:

the angle (v, u), (v, w) does not lie inside the interval (π/2ε, π/2 + ε); experimentally, we determined ε = 0.05 rad; the line u, w lies completely inside the segmented volume.

The second condition detects the paths which lie inside the same tube as the parent path, however, the detection is applied only to paths which are not almost perpendicular to the parent path. The skeleton constructed by the graph-based method presented above cannot be used directly as the domain for the finite element method based on the beam formulation, which requires smooth geometrical representation where every node is equipped with both positional and rotational degrees of freedom (DoF).

Therefore, the next step in the reconstruction process is based on Bézier curve fitting: a set of cubic Bézier curves is fit to each segment of the skeleton using the recursive algorithm presented in [START_REF] Schneider | An algorithm for automatically fitting digitized curves[END_REF]. Then, the 6 DoF nodes are sampled along the Bézier curves adaptively: the density of sampling increases along the segmented with higher curvature in order to improve the quality of the discretization.

Patient-specific composite model

Similarly as in [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF], we consider the registration as a dynamic process. This avoids having to set Dirichlet boundary conditions such that the stiffness matrix is invertible. Such boundary conditions would not make sense in an ICP-based approach. The dynamic system of equations is given by

M ü + B u + K(u) = f (1)
where M is the mass matrix, B is the damping matrix, K is the stiffness matrix, u is the vector of nodal displacements and f is the vector of external forces applied to the deformable object.

In this paper, the liver is modeled using a composite approach which takes into account the mechanical behavior of the three components of the liver: parenchyma, vessels and Glisson's capsule. Since the mass of both the vessel walls and the capsule is negligible w. r. t. the entire body, the mass matrix M is based only on the mass of the parenchyma. In our approach, the damping matrix B is obtained by an approximation

B = r M M + r K K (2)
where r M and r K are Rayleigh mass and damping, respectively. Finally, the mass matrix K is the part of the equation which is composed from three different contributions corresponding to the three components of the liver as explained in the following section. FE formulation of the liver parenchyma While time-dependent phenomena related to viscosity are not considered in this work, we employ co-rotational finite elements to model the parenchyma relying on the geometrical co-rotational method proposed in [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF]. Although the technique relies on linear stress-strain relationship, large displacements including rotations can be considered as correctly approximated. Denoting p a generic P1 tetrahedral element of the parenchyma, the 12×12 element stiffness matrix is then computed as

K p = R p { V p B p D p B p dV }R p (3)
where B p is the strain-displacement matrix, D p is the stressstrain matrix and R p is a rotation matrix derived in the corotational formulation. The details related the matrices in Eq. 3 can be found in [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF]. While both B p and D p are constant during the simulation, R p must be updated in each step of the simulation.

In the case when lesions (such as tumors) are present inside the parenchyma, these can be taken into account also from the mechanical point of view as they usually introduce significant heterogeneity. This heterogeneity can be straightforwardly included in the model since the employed meshing technique [START_REF] Boltcheva | Mesh generation from 3d multi-material images[END_REF] allows for creating heterogeneous meshes according to different labels in the input segmented map.

Then, various constitutive parameters can be assigned to the elements that belong to the area representing the lesion.

FE formulation of the vascular structures

Similarly, the co-rotational approach is used to model the vascular structures. The serially-linked beam elements based on Tymoshenko formulation are employed, taking into account the hollow structure of the vessels via moments of inertia. As the beam formulation considers both positional and rotational degrees of freedom (allowing for modeling the twists and torques), each beam element v is modeled with a 12×12 elements stiffness matrix K v . Although the size of the local system is the same as in the case of parenchyma, the matrix has a completely different structure, as it describes an element given by two nodes each having six degrees of freedom. Thus, it relates angular and spatial positions of each end of a beam element to the forces and torques applied to them. The exact definition of the matrix K v can be found in [START_REF] Duriez | New approaches to catheter navigation for interventional radiology simulation 1[END_REF].

FE formulation of the Glisson's capsule

Due to the thickness of the Glisson's capsule (<20µm), it is not possible to model such a thin structure with classical tetrahedral elements in real-time as the model would require an extremely dense mesh. Instead, modeling the capsule with two-dimensional elements that abstract from the thickness in the third dimension seems to be a natural choice. In the elasticity theory, this functionality is usually provided by membrane elements. To maintain simplicity of the composite model we choose constant strain triangular (CST) elements based on co-rotational formulation. For each CST element c, the 9×9 stiffness matrix K c is given as

K c = V B c D c B c dV (4) = h A B c D c B c dA (5) = hAB c D c B c (6)
where B c is the strain-displacement matrix, D c the material matrix, h is the thickness and A the area of the element. In the previous [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF] follows from the fact that we assume constant thickness of the element and ( 6) follows from the fact that the strain-displacement matrix is constant in our case. The strain-displacement matrix for the CST element can be expressed as:

B c = 1
The values x i j = x ix j and y i j = y iy j are computed from the x or y coordinates of the nodes i, j of the triangular element. The reader can refer to the respective literature [START_REF] Felippa | A study of optimal membrane triangles with drilling freedoms[END_REF] for more thorough description.

Complete liver model

The complete model of the liver takes into account all three components described above by combining the contributions of parenchyma, vessels and capsule. The coupling between the vessel elements (beams) and parenchyma elements (tetrahedra) is described in detail in [START_REF] Peterlik | Modeling and real-time simulation of a vascularized liver tissue[END_REF]. Briefly, each vessel node (having 6 degrees of freedom to account for torques in vessels) is coupled with a tetrahedra in which it is located via barycentric coordinates. This coupling remains constant during the simulation and can be described via matrix J v→p which is the Jacobian matrix of the coupling.

The coupling between the parenchyma and capsule is straightforward, since the triangles used as the domain for the CST formulations are the surface faces of the volume mesh. Therefore, for a given triangle with vertices v 1 , v 2 and v 3 , the corresponding tetrahedron (sharing the same three vertices) is found and a 9×12 permutation matrix P c→p mapping the triangle vertices to tetrahedra vertices is computed.

Without loss of generality, let us suppose a tetrahedral element e which receives a mechanical contribution from both the capsule and vessels. The composite element stiffness matrix K e is then computed as

K e = K p + J v→p K v J v→p + P c→p K c P c→p . (8)
The beam and triangular elements together with the areas of the parenchyma having different Young's modulus introduce heterogeneity and anisotropy into the simulation of the organ. Therefore, either direct solver or preconditioners must be used to solve the system in each step of the simulation.

Numerical solution

The dynamic system given by Eq. 1 is integrated using implicit backward Euler scheme presented in [START_REF] Baraff | Large steps in cloth simulation[END_REF]. The scheme requires to compute the inversion of the system matrix in each time step. Since the system matrix consists of the contributions from the vessels and capsule which introduce heterogeneity and anisotropy, the convergence of iterative solvers such as the conjugate gradient (CG) is jeopardized [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. In our case, due to the heterogeneity and anisotropy introduced by the vessels, the CG requires an important number of iteration or might not converge. Therefore, the system matrix is factorized using the LDL decomposition. In the future, we will use the technique based on the asynchronous preconditioning [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF] in order to improve the computational time.

Image-guided Non-rigid Registration

To register the pre-operative liver model onto the intra-operative laparoscopic video stream, we rely on 1) a reconstruction of a three-dimensional point cloud from the stereoendoscopic camera view, and 2) a multi-step registration process. The main requirements and constraints of this problem can be summarized as follows. First, the stereo-endoscope is capable of capturing only a part of the liver (about 30% to 40% of the whole surface), and we must rely on this partial information to drive the registration process. Second, the liver sustains significant deformations due to the inflation of the abdomen with CO2, and the reference frames from the pre-operative CT scan is completely different from the reference frame of the stereoscopic image, leading to a very poorly initialized registration problem. Third, the liver is continuously moving, under the influence of both breathing and cardiac motions. The following sections describe the approach we propose to address these different points.

Atlas-based transfer of anatomical features

In order to improve the accuracy and robustness of the registration process, we augment the pre-operative liver model with anatomical features. The first step consists of labeling corresponding areas on both the pre-operative mesh and the point cloud to ensure that this point cloud, corresponding to the visible part of the liver, is accurately matched on the pre-operative surface mesh. Three anatomical features used as landmarks: the umbilical notch where the right margin of falciform ligament is attached to the liver, the anterior margin of the liver (see Fig. 3), and the vena cava which is tightly connected to the posterior part of the liver. The second step consists in assigning different roles to these landmarks. Some of them are used to drive the initial registration, others are used as boundary conditions during the realtime registration process. The anterior margin and umbilical notch are used in the registration process, and the vena cava as a boundary condition of the biomechanical liver model. Since the anterior margin of the liver can be easily identified in most human liver, it is automatically detected as follows: the edges separating two triangles with sufficiently different directions of normals are selected as seeds. Then, the ridge line is extended from the seed edges; if no extension is found, the seeds are removed. The definition of the two thresholds, one for the seeds and one for the extensions, is based on statistics on the normals for the whole mesh. Iteratively, the ridge line is reconstructed.

Depending on the image modality, the two other anatomical features are either difficult or simply impossible to segment from the pre-operative data. To solve this issue, these features are transferred onto the three dimensional mesh via a statistical atlas hand-built by experts. Details regarding the construction of the atlas are described in [START_REF] Plantefeve | Atlas-based transfer of boundary conditions for biomechanical simulation[END_REF]. For each liver used to create the atlas the segmentation and the landmarks labelling require anatomical expertise and can take several hours. A set of pre-operative medical images of the organ of interest is segmented manually by an expert to obtain segmented map of the organ and anatomical landmarks. From each segmented maps, surface meshes, one for the whole surface and one for each structure of interest, are generated using the method described in section 2.1.1. First, a number of points is sampled onto the surface of each mesh. The same number of points and numbering is used for each point cloud corresponding to the same structure. Then, all the point clouds are aligned in a common reference frame using a generalized Procrustes analysis (GPA) [START_REF] Gower | Generalised Procrustes Analysis[END_REF]. The GPA determines for each point cloud the similarity transformation SIM i providing an optimal alignments of all the models in the database:

SIM i : S i → C SIM i = s i R i T i (9)
P i → P i (10)
where S i is the space of the i-th segmented model, C the common reference frame, P i the initial point cloud, P i the aligned point cloud, T i the translation which align the center of mass of the ith point cloud and the origin, R i a rotation, and s i a scale factor. Both R i and s i are computed iteratively using the matrix P i P , where P is the mean shape recomputed at each iteration. We have:

SVD(P i P ) = U i Σ i V i (11)
R i = U i V i (12)
s i = tr(P i P ) n||P i || F || P || F (13)
where SVD(.) denotes the singular value decomposition, tr(.) the trace of a matrix, and ||.|| F the Frobenius norm of a matrix. Then, the principal component analysis (PCA) is performed on each structure of interest to compute the principal modes of deformation across the database. We keep only the most significant modes, i. e. the modes responsible for more than 90% of deformations. Assuming the Gaussian distribution of the modes, we can determine the mean position and the standard deviation associated with each structure.

The atlas is finally registered on a three-dimensional patientspecific liver model using the physically-based registration method described in section 2.2.3.

Intra-operative 3D reconstruction

The estimation of the three-dimensional shape of the liver from intra-operative images can be performed in different ways depending on the image modality [START_REF] Maier-Hein | Optical techniques for 3d surface reconstruction in computer-assisted laparoscopic surgery[END_REF]. In this work, we make no assumption about specific acquisition technologies available in the operating room (such as intra-operative MR or CT scanner, or depth camera), but we rely only on images provided by the stereoscopic camera for the estimation of the three-dimensional shape of the liver. Although this technology is not currently widespread in operating theaters its usage is rapidly growing. The only requirement is that the camera be calibrated at the beginning the intervention which takes only few seconds.

Let us assume the stereoscopic pair of images I l and I r . We extract from this stereo pair points of interest (i. e. features) that are sufficiently reliable for the 3D reconstruction. Several detectors have been reported in the literature (e.g. in [START_REF] Gauglitz | Evaluation of interest point detectors and feature descriptors for visual tracking[END_REF]), we employ the Speeded-up Robust Features (SURF); our choice is motivated by its particular suitability for robotic-guided endoscopy applications [START_REF] Elhawary | Robust feature tracking on the beating heart for a robotic-guided endoscope[END_REF].

When using this detector on the image pair (I l , I r ) with an appropriate threshold, we obtain two sets of features E l = (x 1 l , . . . , x n l ) and E r = (x 1 r , . . . , x m r ) where it is necessary to estimate for each feature

x i = (u i , v i ) the 3D point X i = (U i ,V i ,W i )
. This is done by establishing a correspondence between image points x l ←→ x r , using a descriptor-based matching method [START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF]. Once a correspondence is found, a sparse set of m 3D points, denoted T = (X 1 , . . . , X m ), is reconstructed using triangulation algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] and a surface mesh S is interpolated using Moving Least-Square surface approximation [5,?].

The point cloud has to be labelled with the anatomical landmarks described in section 2.2.1. An operator selects (e. g. by clicking) on one of the stereoscopic images some points on the border of the three regions of interest, the umbilical notch, the anterior margin, and the liver surface, to form closed areas. This selection takes only few seconds. Once the targeted regions are selected, we build for each of them a three-dimensional point cloud using the method described above. These point clouds are labelled according to their corresponding regions as illustrated in Fig. 3a.

Elastic registration

At the beginning of the registration process, the biomechanical model (source) is in its pre-operative configuration and needs to be deformed to match the intra-operative point cloud (target). This match is only partial since the target represents only about 30% to 40% of the total surface of the liver. The biomechanical model described in section 2.1.2 is used to simulate the behaviour of the tissue while the registration constraints imposed to the deformable object are modelled with penalty forces f ext . Theses penalty forces can be seen as a mean of imposing a displacement with a tolerance. Therefore, the result of the registration does not depend on the absolute values of the stress parameters but rather on their relative ratio, as already shown in [START_REF] Wittek | On the unimportance of constitutive models in computing brain deformation for image-guided surgery[END_REF], and on the relative intensity of external forces versus internal forces. The goal of the registration process is to minimize the energy function associate to theses forces. At each time step, f ext is recomputed as follows. First, we define a set of corresponding pairs {p i , q i } where p i is a point of the point cloud and q i is a point of the biomechanical model surface. To compute these pairs, we project each point of the point cloud with label j onto its closest neighbour on the part of the surface mesh with the same label using the method described in [START_REF] Plantefeve | Automatic of pre and intraoperative data using anatomical landmarks for augmented laparoscopic liver surgery[END_REF]. Then, for each pair we define an external force:

f ext i = k(1 -arctan(s|| --→ p i q i ||)) --→ p i q i || --→ p i q i || ( 14
)
where k is a scalar stiffness coefficient (in [N/m]) and the term 1arctan(s|| --→ p i q i ||) is an asymptotic penalty function of the distance (in [m]) which includes also the scale factor s. The stiffness coefficient k is defined in three different ways depending on the type of the corresponding feature: thus, k is one of k un , k am and k surf where the identifiers stand for the umbilical notch, anterior margin and surface, respectively. At the beginning of the registration, the reliability of matching associated to each type of feature is reflected in the corresponding values of the stiffness coefficients: k un > k am > k surf . As the registration proceeds, k un remains constant, whereas both k am and k surf evolve in order to prevent the energy minimization process from falling into a local minimum. Therefore, during the registration, it is detected when the mechanical system has reached a plateau, and the stiffness coefficients k j for j ∈ {am, surf} are updated as

k j min + 2(k j max -k j min )( n r j ) 2 ( n r j ) 2 + 1
for n < r j and k j max for n ≥ r j [START_REF] Plantefeve | Automatic of pre and intraoperative data using anatomical landmarks for augmented laparoscopic liver surgery[END_REF] where n ≥ 0 is the plateau index (which is incremented each time a new plateau is detected) and r j > 0 an integer which controls the increase rate of k j . The registration is stopped as soon as n = max(r j ) + 1.

After attaining convergence, the temporal registration step is performed. The temporal registration permits to register the biomechanical model at each frame acquired from the endoscopic camera, following the motion of the liver (which moves due to the heart beating and/or respiratory motion). As a good initialization of the three dimensional model configuration is available, the point cloud is down-sampled and the labels (umbilical notch and anterior margin) are no longer used. The motion of the liver is captured using an optical flow algorithm [START_REF] Bouguet | Pyramidal Implementation of the Lucas Kanade Feature Tracker: Description of the algorithm[END_REF] that tracks image-points from the liver surface based on a brightness consistency constraint 4. The combination of optical flow and SURF features has proven to be robust to track heart motion in laparoscopic images, where real-time is needed [START_REF] Elhawary | Robust feature tracking on the beating heart for a robotic-guided endoscope[END_REF] and was successfully translated to liver motion estimation [START_REF] Haouchine | Towards an accurate tracking of liver tumors for augmented reality in robotic assisted surgery[END_REF]. Indeed, at this stage, a real-time simulation is mandatory since video frames are acquired at least at 30 Hz and the registration w. r. t. the nth frame must to be performed before the frame n + 1 is available. Furthermore, if the boundary conditions were not used during the initial registration because of the rigid transformation between the pre and intra-operative liver model, we activate them during the temporal registration. However, some of the liver boundary conditions are mobile, such as the ligaments that attach the liver to the diaphragm. For this reason, we only use the vena cava and the entrance points of the liver vascular tree as boundary conditions to allow the model to move with the diaphragm.

Results

Throughout the various stages of the pipeline, errors arising from different sources may degrade the quality of the final result. First, the segmentation of the liver and its vascular network is operator-dependant and sensitive to the quality of the medical images. The error arising from the segmentation is not easy to evaluate but a lower bound is given by the voxel size which in our case is 0.85 mm x 0.85 mm x 0.70 mm. Second, the quality of the three-dimensional mesh generated from the segmented maps depends notably on the number of its vertices. We generated a high-quality mesh and computed the Hausdorff distance between this mesh and the mesh that we use for computation. We found a mean Hausdorff distance of 0.2 mm and a maximal error of 4.9 mm. Third, the point cloud quality is subject to errors arising from the calibration of the camera, the image quality (sharpness, texture. . . ) and from the triangulation method (features detection, matching. . . ). As we do not have any ground truth for the laparoscopic data, it is really difficult to estimate the error in millimetres between the point cloud and the real position of the organ. Fourth, the use of an atlas to transfer some landmarks onto the three dimensional model may lead to some error in their location. We tried to evaluate this error and its influence in section 3.1. Finally, we evaluate the quality of the registration using our complete biomechanical model of the liver in section 3.2.

Transfer of Boundary Condition

In this section we assess the quality of the boundary condition transfer and its influence on the biomechanical model behaviour. The atlas was created using ten liver models and their respective anatomical features obtained from segmented abdominal CT scans. To ensure the quality of the atlas, we evaluated the standard deviation of the umbilical notch and vena cava positions over the samples; we obtained σ un = 11.12 mm and σ vc =16.56 mm respectively. The average variability represents less than 9.7% of the size of the organ showing a strong consistency among the feature positions. We registered the atlas onto several livers: ten from the atlas following the left-one-out principle and two livers not included in the atlas (sample #11 and #12). The registration errors are presented in Fig. 5.

The results show that the surface registration is very accurate and all the boundary condition location errors are under 10 mm except for one sample. This value is less than the standard deviations in the atlas. To evaluate the influence of errors in boundary condition location, we altered the position of the boundary conditions compared to the reference: the falciform ligament attached to the umbilical notch and the vena cava. We used a liver segmented in the supine position to set up a simulation where gravity is re-oriented in order to simulate a deformation of the organ in the flank position. This simulation was performed using both the reference boundary conditions location and the altered ones. Then, we compared the final shape of the liver simulated with altered boundary condition to the final shape of the reference liver. The results are presented in Fig. 6(a). They show that small displacement (below 3 mm) lead to an error below 1 mm, however, the relation between the displacement and the error seems to be quadratic. Nevertheless, a deeper analysis is needed to confirm this tendency. The same kind of experiments was conducted to evaluate the influence of errors in the boundary condition stiffness, the results are shown in Fig. 6(b). The values of the reference stiffness were set according to the ligament and vena cava Young's modulus found in the literature [START_REF] Umale | Characterization and modeling of abdominal organs[END_REF] (20 MPa for both) and their geometry. The results show that the influence of the stiffness is very small when compared to the influence of the location of the boundary conditions. This can be justified since the ratio between the liver internal forces and the external forces imposed by the boundary condition is very low. Therefore using twice or half the reference stiffness does not change this ratio significantly.

Patient-specific modeling and registration

The aim of the method is to estimate the location of internal structures of the organ which undergoes important intra-operative deformations when compared to its initial pre-operative configuration. Validating the registration in this context is very challenging, since the optimal ground truth would be a 3D reconstruction of the organ at the intra-operative stage, which requires the intra-operative CT or MR scan of the patient. Access to these techniques is very limited and almost impossible to use on human subjects. Consequently, we decompose the validation process in two parts: first, we present quantitative results on in silico data (we decided not to use a phantom for an intermediate validation because of the low complexity of such model and non-realistic conditions for surface patch reconstruction). Second, we report qualitative results using real endoscopic data.

Patient-specific modeling

A heterogeneous anisotropic biomechanical model was generated from the patient image data using the method described in section 2.1. In order to obtain the ground truth for the registration and to measure the influence of the biomechanical model used for the deformation, we applied a pressure on the model surface to simulate the effect of the pneumoperitoneum. We used two parameter sets; one representing a healthy liver (parameter values are set as reported in [START_REF] Umale | Experimental in vitro mechanical characterization of porcine Glisson's capsule and hepatic veins[END_REF]):

E parenchyma = 3.5 kPa, ν parenchyma = 0.45, E vessels = 620 kPa, ν vessels = 0.4, E capsule = 8000 kPa, ν capsule = 0.45, and the other a cirrhotic liver: E parenchyma =30 kPa. All the other parameters are identical for both cases. We also deformed two homogeneous models of the parenchyma using the Young's modulus and the Poisson's ratio of the healthy and cirrhotic liver. Then, we measured the distance between the final states obtained after deformation. We compared the healthy liver to the cirrhotic one and the homogeneous models to their heterogeneous counterparts. Results presented in Fig. 7 show that the vessels and the capsule influence the deformations, particularly in the case of the healthy liver where the mean Target Registration Error (TRE) between the final configuration of the heterogeneous and homogeneous model is 10.1 mm. This error represent 43% of the homogeneous healthy liver deformation. The same analysis performed for the cirrhotic liver quantifies the influence of the vessels being 18% w. r. t. the entire deformation.

Registration

As only a part of the liver surface is visible to the endoscopic camera during a surgical intervention, we extracted a portion of the deformed surface from the deformed configuration, obtained with the healthy and cirrhotic heterogeneous liver model. We kept two deformed configurations having an area ranging from 50% to 10% of the whole surface and evaluated the registration robustness w. r. t. the amount of visible surface.

The undeformed healthy and cirrhotic biomechanical models were subsequently registered onto these partial surfaces. Then, we compared the shape of the models after registration with the deformed configurations of the whole meshes. We used a tetrahedral mesh composed of 2193 tetrahedra and 547 vertices. The Rayleigh mass and the Rayleigh damping coefficients were set to 0.1. For all the registrations we adjusted the parameters according to the Young's modulus value. In this manner, we ensured that the ratio between the external and the internal forces remains the same independently of the actual value of the Young's modulus. We use a factor 10 -3 for k un , k am max , and k surf max , a factor 10 -4 for k am min and we set k surf min = 0, r am = 2, and r surf = 10. The computational time needed for the initial registration remains below 30 s which is acceptable for the given scenario.

We computed the TRE on the mesh surface and in the mesh volume to measure the accuracy of the registration. The results are presented in Fig. 8. We obtained a mean error below 3 mm in all cases except when the visible surface represent less than 20% of the liver surface, which would correspond to a very limited surface reconstruction.

We also tested the influence of our heterogeneous liver model versus a simplified homogeneous liver model on the registration. In the case of a cirrhotic liver the registrations of the simplified model leads to similar TRE showing that the vessels have a low influence on the model. However, in the case of a healthy liver the volume TRE increase when the registration is performed with a simplified biomechanical model. Additionally, we studied the influence of the point cloud size on the registration results. Provided that the density of points remains above 0.7 points per square centimetre, the mean TRE between the registration result using the complete point cloud (3.9 points per square centimetre) and the down sampled one does not exceed 1 mm. Moreover, as the umbilical notch position is transferred to the patient liver via an atlas, we conducted a preliminary study to evaluate the influence of the umbilical notch position on the registration results. We moved the umbilical notch around its real position on the liver surface. The TRE between the registered source mesh using the reference umbilical notch position versus the altered position highly depends on the displacement direction. The worst scenarios occurs when the umbilical notch is displaced along the anterior margin towards the right lobe. In this case the mean TRE is increased by 3.5 mm for a position shift of 17.5 mm (mean hausdorf distance of 4.8 mm). However, for all tested scenarios, the TRE remains below 2 mm as long as the displacement of the ligament does not exceed 10 mm (mean hausdorf distance of 1.0 mm). Nevertheless, additional experiments would be necessary to conduct a proper analysis of the impact of position error of the umbilical notch on the registration result.

Finally, we tested our method on in-vivo human data using two sets of patient images. As in-vivo a direct quantitative evaluation is not possible, we perform only a visual qualitative assessment presented in Fig. 9. Nonetheless, even in this case, it is challenging to provide comprehensible vi-sualization of the registration results. However, the point cloud was accurately matched to its corresponding part onto the three dimensional mesh and the mean Hausdorff distance between the surface mesh and the point cloud was below 1.1 mm.

Conclusion

We presented a complete pipeline to register a pre-operative patient-specific three-dimensional model of the liver to sparse and incomplete intra-operative data extracted from stereo-endoscopic images. This pipeline includes a semiautomatic segmentation of the liver and its vascular system, an automatic generation of the biomechanical model from the segmented maps, a real-time simulation of the complete liver model. It also involves the automatic adaptation of boundary conditions using a statistical atlas, and an ICPbased non-rigid registration process. We performed a quantitative assessment of the registration accuracy using the synthetic data for which the ground truth is known. The evaluation was done for two different sets of parameters corresponding to healthy and cirrhotic liver. In the case of the healthy liver, the influence of the vascular network on the overall deformable behavior was more significant than in the case of the cirrhotic liver which is much stiffer. In both cases, the elastic registration was capable of providing a good estimation of the actual position of the internal structures from the surface data: we showed that only 30% of the organ surface captured by the camera is sufficient to reduce the mean TRE measured both on the surface and in the volume below 3 mm. This value is smaller than the current laparoscopic procedure safety margins which are usually between 10 mm and 25 mm around the possible location of a tumor. We have also performed a qualitative assessment using real patient data. Both the evaluation using the synthetic and the real data demonstrated that using a biomechanical model results in good accuracy and robustness of the method. We are currently developing a prototype of our system that can be deployed in the operating room. This requires to further automate the pipeline described in this work, and to reduce the number of parameters involved in the process. We are also working on a more thorough validation protocol based on an intra-operative CT scan in which the position of the camera could be determined as well as the location of the tumor(s). This will allow for a direct comparison of the estimated and actual target locations.

In the future, we plan to improve the definitions of the boundary conditions for the temporal registration phase which is an important and challenging problem. We will also try to design an in vivo validation setup to obtain quantitative errors on real data. In addition, we plan to improve the visualization of the vascular network superimposed on the laparoscopic view to ensure better depth perception.

Fig. 2 :

 2 Fig. 2: Scheme showing the pipeline used for the model reconstruction: (a) input CT or MRI data, (b) segmented vascularization, (c) segmented parenchyma, (d) skeletonization of vascular trees, (e) meshing of segmented volume and (f) composite model composed of the parenchyma mesh and beam elements sampled from Bézier curves fitted to the vascular skeleton.

2A   y 23 0 y 31 0 y 12 0 0 x 32 0 x 13 0 x 21 x 32 y

 32 23 x 13 y 31 x 21 y 12  

Fig. 3 :

 3 Fig. 3: Anatomical landmarks. On the laparoscopic image: the liver surface (in green), the anterior margin (in orange) and the umbilical notch (in yellow).

Fig. 4 :

 4 Fig. 4: A SURF feature being tracked over time (600 frames) using LK optical flow. (top) tracking performed in both a forward and a backward time direction. Graph shows X and Y coordinates as a function of video frames. (middle) The selected feature located in the endoscopic images at selected frames. (bottom) 50 × 50 pixel window centred around the feature location at every 100 frames.

Fig. 5 :

 5 Fig.5: Mean Hausdorff distances between the segmented liver surface and the registered mean shape (in red) between the manually segmented umbilical notch and the umbilical notch transferred from the atlas (in green) and between the manually segmented vena cava and the vena cava transferred from the atlas (in blue).

Fig. 6 :

 6 Fig. 6: Sensitivity study. Mean Hausdorff distance between the final configurations of a reference model and a biomechanical model with altered boundary conditions position (a) and stiffness (b).

Fig. 7 :

 7 Fig. 7: Surface and volume TRE between deformed configurations. A) heterogeneous liver model: healthy vs. cirrhotic , B) healthy liver: homogeneous vs. heterogeneous model, C) cirrhotic liver: homogeneous vs. heterogeneous model.

Fig. 8 :

 8 Fig.8: Surface error distributions w. r. t. the percentage of visible surface. The boxes with solid lines represent the results for the healthy liver and the dashed boxes correspond to the cirrhotic liver (the bar-and-whisker graph shows the mean, upper and lower quartiles, and the maximal and minimal error). For the partial surface representing 10% of the liver surface we also measured the errors on a sub-part corresponding only to the volume under the visible surface from the camera viewpoint (10p).

Fig. 9 :

 9 Fig.9: Non-rigid registration between intra-operative and pre-operative data using in-vivo human data. The overlay of the liver surface and the vascular network permits to guide surgeon during the operation.


D Ferraro

Y Lin

B Teste

Damien Talbot

L Malaquin

S Descroix

Ali Abou

A Abou-Hassan

Continuous chemical operations and modifications on magnetic γ-Fe2O3 nanoparticles confined in nanoliter droplets for the assembly of fluorescent and magnetic SiO2@γ-Fe2O3

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Due to their specific properties rising from their small dimensions, nanoparticles (NPs) are used in various applications such as catalysis, optics, electronics and biomedical diagnosis or nanomedicine. 1,2 To extend even more the range of properties and applications of NPs, several nanomaterials can be combined into one single nanostructure to form a multifunctional nanohybrid. This can be accomplished by controlling the assembly of smaller nanostructures (building blocks) of various shapes and compositions to create complex and single nanoscale object, having different properties provided by each of the different nanomaterials.

3

In the context of nanoparticle synthesis and assembly, different microfluidic strategies had already been shown with a high level of synthesis control enabling the production of high quality nanostructures with single or multi-functionalities. [START_REF] Song | [END_REF]5 Particularly, in materials science, droplet microfluidics has been identified as a powerful tool to control chemical reactions through their compartmentalisation inside individual droplets, where each droplet acts as a micro-reactor allowing reaction parallelisation, reduction of reagents volume and accurate mixing associated to an enhancement of mass and heat transfer compared with bulk reactions. 6 Moreover, it is also possible to manipulate individually each droplet to perform complex operations such as splitting, merging, sorting that allow to combine different chemical reactions. 7,8 Despite several convincing reports on the use of microfluidics for the elaboration of nanomaterials through chemical synthesis or self-assembly, [9][10][11] there is still rooms for high impact studies relying on technology for complex workflow implementation at microfluidic scale. In particular, manipulations of magnetic nanoparticles in microfluidic systems are still limited to single step reactions, 12 or multistep continuous modifications, 13 when operations as extraction, washing or NP functionalisation are not required or performed out from the microfluidic devices by conventional approaches. This is mainly due to their nanoscale dimension which renders difficult their handling in confined condition, where low amount of NPs, or in other words low amount of paramagnetic materials, results in low magnetic force. In these states typical macroscale protocols cannot be directly exploited.

Herein we provide for the first time an automated microfluidic platform capable of manipulating nanoliter droplets and performing multistep and automated sequences of operations on magnetic NPs. These operations include screening of colloidal interactions magnetic decantation, colloidal redispersion, washing, surface functionalization, heating and colloidal assembly. In order to illustrate the potential of this microfluidic technology we report the elaboration of fluorescent and magnetic nanoassemblies from fluorescent silica and magnetic NPs as building blocks. Nanohybrid structures with dual functionalities are very promising and have proven their potential in different biomedical contexts such as drug delivery, magnetic resonance imaging (MRI) and hyperthermia cancer treatment. , offers the possibility to generate a deterministic and confined train of droplets in a range between 50 and 300nL (see example in Figure 1b) with an error in volume of about 2%. In order to generate stable Please do not adjust margins Please do not adjust margins droplets, FC-40 is mixed with 2% (w/w) of a fluorinated surfactant (1H,1H,2H,2H-perfluorodecan-1-ol, by Fluorochem). Additionally, the microtiter plate is entirely covered with a small layer (about 2mm) of mineral light oil (Sigma Aldrich) to avoid the reagent evaporation. Since the syringe pump can work in both injections and withdrawal modes, droplets can be moved forward and backward in the capillary. This capability will be exploited for inducing a recirculation inside the droplet, which helps the internal mixing. 16 Once the train of droplets is generated, droplets can be driven between the magnetic tweezers to achieve the manipulation of the magnetic NPs and through a homemade heater (see ESI, Figure S1), which allows the temperature control of the droplet. The magnetic tweezers consist of two sharp magnetic tips (AFK502 Imphy Alloys) facing each other across a capillary where the droplets are flowing. Each tip of the tweezers can be activated by a magnetic coil. Therefore it is possible to switch ON and OFF the magnetic field simply by controlling the coil current. The power of this technology had already been validated for extraction or agglutination of magnetic micrometric beads for protein detection, 17,18 however the manipulation of colloidal magnetic NPs has never been demonstrated before in droplet microfluidics.

The challenge of the manipulation of NPs compared to micrometric beads in microfluidic droplets could be explained considering the forces acting in the system. For example, it has been shown that to extract micrometric particles from a droplet, the magnetic field allows to collect them in a cluster that can be extracted if the magnetic force is higher than the capillary force 17 (see Movie 1, ESI). After that, the cluster of beads can be merged and released in the following droplet. However, in case of ferrofluid the small dimension of the particles favours the dispersion forces relative to Brownian motion. In particular, in this case we noticed that the entire aqueous droplet is trapped, but no cluster formation and extraction are observed (see Movie 2, ESI). In fact, NPs experience two opposite forces: the magnetic force which tends to aggregate them and the colloidal force which keeps the suspension stable, resulting in a metastable system where the droplet is trapped at the magnetic tweezers (see paragraph 2 in ESI). Nevertheless, it is interesting to notice that even if it does not permit particle extraction, this effect can be advantageously exploited for merging two following droplets. Size sorted cationic ferrofluids of maghemite (γ-Fe 2 O 3 ) NPs suspended in nitric acid (pH = 2), with different diameters (d 0 between 5-9nm; see ESI, Figure S3 for more details) were synthesised according to our previous described methods. 19,20 A phase diagram describing the efficiency of ferrofluid droplet trapping effect was built to report the influence of NPs diameter and iron concentration (Fe) at constant oil flow rate (Q=0.02L•s -1

) (see ESI, Figure S4). Q was chosen in order to be able to trap also the droplet of ferrofluid with smaller NPs at lower concentration of iron. Since the magnetic force is proportional to the quantity of magnetic material, as expected bigger and higher concentrated NPs are more suitable for the droplet trapping. For example fixing Fe=0.

3 mol•L -1
, only droplet containing NPs having d 0 ≥7nm can be trapped. In order to perform complex ferrofluid manipulation, it is necessary to be able to extract the NPs from the ferrofluid. This can be achieved increasing the magnetic field to overcome the colloidal force and generate a cluster, as described for micrometric particles, or destabilizing the colloidal suspension. Considering the difficulties to generate high magnetic fields using electromagnets, the latter solution was found to be simpler and more straightforward to implement experimentally. In particular, this was achieved by merging the droplet of ferrofluid with a droplet of acetone, which increases the solvent permittivity leading to NPs flocculation. At this point, NPs are not stable in the solutions, they aggregate quickly and accumulate at the tweezers forming a cluster that can be extracted from the initial droplet (see Movie 3, ESI) and released in any following droplet or washed with subsequent droplets.

Based on these validations, we demonstrated the potential of the platform to perform series of complex synthesis workflows. In particular we focused on the elaboration of magnetic and fluorescent nanoassemblies from maghemite (γ-Fe 2 O 3 ) and silica NPs as building blocks. We have selected γ-Fe 2 O 3 NPs of d 0 90,23 nm as the magnetic building block since NPs with this size are easier to manipulate using our tweezers and have proven their efficiency in different applications related to nanomedicine. 21 The other building block is amine functionalised fluorescent silica NPs (RITC-SiO 2 ) (d2004 nm; ESI Figure S5) that have been synthesised according to the modified Stöber method using the Rhodamine B isothiocyanate (RITC,  absorption = 555 nm) as a fluorescent dye. 22 The amine functionalised silica NPs are stable at pH= 7 in 10 To build our assembly through electrostatic interaction at pH=7 between the positively charged silica and the maghemite NPs, we Please do not adjust margins Please do not adjust margins decided to coat their surface with citrate molecules since maghemite NPs are not stable at this pH. 13 The role of this coating is to ensure their colloidal stability due to carboxylate functional groups on their surface, which also guarantee a strong interaction with positively charged silica NPs. Therefore, the workflow necessary for reaching the final assembly, which is completely performed in droplets, is illustrated in Figure 2a. At first, the train of droplets sketched in Figure 1b is ). Based on the trapping mechanism previously described, the droplet containing -Fe 2 O 3 is trapped at the tweezers and merged with the following droplet containing trisodium citrate (step 1, Figure 2b-g) to perform the NPs citration. The addition of trisodium citrate solution induced nanoparticle destabilisation and precipitation in the droplet (see Movie 4 ESI) due to the high ionic strength. To accelerate diffusion and grafting of citrate molecules on the surface of magnetic NPs, the resulting mixture in the droplet was flown in the microfluidic heater introduced above and heated at 60°C for 5 minutes (step 2 in Figure 2). After this step, the resulting suspension of γ-Fe 2 O 3 became optically homogeneous indicating uniform dispersion of stabilized citrate functionalised magnetic NPs in the droplet. Since the sodium citrate is in excess, we expected that the reaction efficiency is 100% and all the NPs resulted functionalized. To remove the excess of citrate molecules, the droplet from step 2 was magnetically trapped again and merged with a droplet containing the mixture of acetone/water, while keeping the magnetic field and the oil flow ON (step 3, Figure 2h-m, Movie 3 ESI). As previously introduced, thanks to the acetone, we could extract the citrated magnetic NPs from the excess of citrate presented in the droplet solution. In the step 4 (figure 2n-s), the NPs cluster was washed with two water droplets, and finally released in the MOPS solution. During washing, the flocculated NPs through which the water droplets flowed were kept trapped at the tweezers. Figure 2n-p and Movie 5 in ESI show how the cluster becomes larger when the water droplets pass through, evidencing slight redispersion of the NPs. Additionally some γ-Fe 2 O 3 NPs are lost in the washing droplet; evaluating the grey levels of the Movie 5, we estimate that this quantity is about the 10% of the total cluster of NPs. Then, when the MOPS solution was mixed with the magnetic NPs (Figure 2q-s and Movie 6 ESI), the resulting droplet was released from the tweezers and flown for 1 minute in the capillary to obtain a better mixing. At this stage, no NPs sedimentation was observed even after 1 hour suggesting the efficient stability of γ-Fe 2 O 3 -citrated NPs in MOPS. Finally, the droplet containing the citrated magnetic NPs dispersed in MOPS was trapped at the tweezers and merged with a droplet loaded with the RITC-SiO 2 fluorescent NPs in order to induce the assembly of RITC-SiO 2 @γ-Fe 2 O 3 nanohybrids through electrostatic interactions. As presented in the previous case, mixing was done flowing the resulting droplet in the capillary for 1 minute in order to benefit from the internal recirculation flows.

Optical fluorescence microscopy was used to confirm the successful formation of RITC-SiO 2 @γ-Fe 2 O 3 assemblies in the nanodroplets. As shown in Figure 3a when the magnetic field (B) is OFF, the nanoassemblies obtained in the droplet are randomly distributed due to Brownian motion; while in presence of a permanent magnetic field (Figure 3b), the fluorescent homogenized nanoparticles lined up (white arrow) along the magnetic field direction (see also Movie 7, ESI) to form chainlike structures due to the attractive interaction between the magnetic moments. Additionally, applying a magnetic field gradient to the RITC-SiO 2 @γ-Fe 2 O 3 assemblies confined in droplet by a permanent magnet, we observed that they tend to agglomerate close to it (see Movie 8, In summary, herein and for the first time, we presented a microfluidic platform for fast (all the process takes about 10 minutes) elaboration of magnetic and fluorescent nanohybrid assemblies based on manipulation of nanoliter microfluidic droplets loaded with nanoparticles and chemicals. This was achieved by using an automated microfluidic platform coupled with magnetic tweezers that allows sequential manipulation and different workflows on ferrofluid suspensions without leaving the microchannel. Our method is versatile; it can be applied to the surface coating of magnetic based NPs with any ligand and moreover to the assembly of different nanomaterials with magnetic NPs. In fact, here microfluidics is not exploited as a production technology, but as a synthesis investigation tool. We believe that magnetic tweezers technology combined with typical advantages of microfluidics could represent a very useful tool for screening chemical conditions saving time and reagents compared with standard chemical approach. However, to reach such a final goal it will be necessary to miniaturise and adapt different techniques known in colloidal chemistry such as light scattering methods in order to acquire fast data on the process and fast screening.

15

 15

  Morpholino)propanesulfonic acid (MOPS) pH buffer and their zeta potential is ζ=+35mV.

Figure 1 .

 1 Figure 1. a) Microfluidic setup for the generation of deterministic trains of droplets. Syringe pump synchronised with a pipetting robot generates the droplet from a microtiter plate, sucking water phase solutions and oil alternatively. These droplets are then driven between the magnetic tweezers (picture in the inset) for ferrofluid manipulation and the homemade heater that can be easily tuned at the desired temperature. b) Train of droplets for nanoassembling production: 1) 50nL of -Fe2O3 solution, 2) 50nL of sodium citrate solution, 3) 200nL of acetone/water solution, 4) 2 drops of 200nL of milliQ water, 5) 50nL of 3-(N-Morpholino)propanesulfonic acid (MOPS) and 6) 200nL of silica fluorescent NPs.

generated: 1 ) 1 ; 2 )

 112 50nL droplet of -Fe 2 O 3 solution (d 0 9nm, γ-Fe 2 O 3 8.2x10 18 NPs•L -1 or Fe 0.44 mol•L -50nL of 0,027mol•L -1 trisodium citrate solution giving a molar ratio Fe/citrate16; 3) 200nL of an acetone/water mixture (50/50 v/v); 4) two water droplets of 200nL each; 5) 50nL of10 -3 mol/L MOPS pH buffer solution (pH = 7) and 6) 200nL of silica fluorescent NPs (d200 nm, RITC-SiO 2  1.88 x 10 15 NPs•L -1

Figure 2 .

 2 Figure 2. (Top) a) scheme of the sequential reaction and workflows for the elaboration of nanoassemblies starting from the aqueous suspension of ferrofluid at pH=2; (Bottom) Manipulation of γ-Fe2O3 NPs confined in droplet by magnetic tweezers: b)-g) (Step 1) γ-Fe2O3 droplet is trapped and mixed with a following droplet which in this case contains sodium citrate solution; (Step 2) the resulting droplet is heated for 5 minutes at 60°C; h)m) (Step 3) destabilizing by acetone/water droplet of the γ-Fe2O3solution, NPs flocculate and can be extracted; n)-s) (Step 4) 2 water droplets pass through the NPs cluster for washing (in this figure we present only the second water droplet) and thenγ-Fe2O3 NPs are released in the MOPS droplet. Transparent droplets are highlighted by red lines on the meniscus. The direction of the flow is from the left to the right.

Figure 3 .

 3 Figure 3. Fluorescent images of the silica nanoparticles coated with iron oxide nanoparticles in the absence (a) and in the presence of the magnetic field (b), oriented as indicated by the white arrow; c-d) TEM images of RITC-SiO2@γ-Fe2O3 nanoassemblies.

Notes and references

We thank Prof. J.C Bacri from MSC lab, Paris 7 for discussions. We thank Dr. Daniel Lévy and Aurélie Di Cicco for their help in TEM measurements. This work was supported in part by the Digidiag project (ANR) from the French government, the ARC foundation.