Stéphanie Marchesseau

Simon Chatelin

Hervé Delingette

Non linear Biomechanical model of the Liver

Understanding and modeling the liver biomechanics represent a significant challenge due to its complex nature.

While many studies have been performed to fit hyperelastic constitutive laws on rheological experiments, they tend to agree about the importance of strain rate in the liver mechanical behavior. Furthermore, as the liver is heavily perfused with blood, its constitutive behavior is greatly porous. Supported by these observations, we developed a porous visco-hyperelastic model as a liver parenchyma material. More precisely, visco-hyperelasticity is obtained through Prony series while the mechanical effect of liver perfusion is represented with a linear Darcy's law. Since this mechanical model is developed in the context of real time surgery simulation, a compromise between biomechanical accuracy and computational efficiency must be found. We propose the Multiplicative Jacobian Energy Decomposition method (MJED) to obtain a fast assembly of stiffness matrices on linear tetrahedral elements. Finally, the relative effects of the hyperelastic, viscous and porous components on the proposed liver model are discussed and compared to some rheological experiments.

Introduction: Clinical context

The liver is among the largest organ of the human body with an approximate weight of around 1.5 kg for adults. The biomechanical behavior of this massive soft organ has been studied on the one hand in the context of car injury [START_REF] Untaroiu | Modeling the biomechanical and injury response of human liver parenchyma under tensile loading[END_REF]. Indeed, liver injuries such as capsule laceration and parenchyma damage are frequently caused by both frontal and side collisions and are associated with high morbidity and mortality rates. The numerous related rheological and simulation studies are therefore performed in a regime of large deformations, high strain rates and are focused on the mechanisms of tissue failure.

On the other hand, the biomechanics of the liver has also been studied in the context of small strain rates for medical applications such as computer-aided diagnosis, therapy guidance, therapy training and therapy planning. For instance, in the context of therapy training, several existing surgery simulators [START_REF] Nickel | Virtual reality training versus blended learning of laparoscopic cholecystectomy: a randomized controlled trial with laparoscopic novices[END_REF] aim at training young surgeons to perform the resection of the gall bladder (cholecystectomy) using minimally invasive surgery. Since the gall bladder sits beneath the right lobe of the liver, those simulators often include a simplified real-time model [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF] of the liver mechanics. There exists several possible therapies that can be performed on the liver related to the presence of hepatic metastases (often related to colorectal cancer) or primary tumors such as hepatocellular carcinoma. Among those procedures, the partial resection of the liver aims at removing the functional regions of the liver that include some tumor lesions. Prior to this surgery, liver biopsies are often performed to examine the nature of the tissue at risk. In both cases, the physician must face the issue of properly localizing the lesions based on pre-operative imaging (CT or MR images) and intra-operative images (often ultrasound images). Additional guidance can be provided to the surgeon by registering the liver shape from its pre-operative configuration to its intra-operative one based on a biomechanical model. In such cases [START_REF] Oktay | Biomechanically driven registration of preto intra-operative 3d images for laparoscopic surgery[END_REF][START_REF] Simpson | Soft Tissue Biomechanical Modeling for Computer Assisted Surgery[END_REF], the mechanical model must cope with large displacement (finite strain) and acts as a data regularizer to constrain the space of deformations. Such a bio-inspired image registration approach was further used to plan the trajectory of biopsy needle in the liver [START_REF] Kobayashi | Soft Tissue Biomechanical Modeling for Computer Assisted Surgery[END_REF], to perform augmented reality to visualize tumor locations during laparoscopic surgery [START_REF] Plantefève | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF][START_REF] Haouchine | Impact of soft tissue heterogeneity on augmented reality for liver surgery[END_REF], to simulate the injection of gas inside the abdominal cavity (pneumoperitoneum) [START_REF] Bano | Simulation of pneumoperitoneum for laparoscopic surgery planning[END_REF].

Anatomy of the Liver

The liver is the heaviest internal organ of the human body, whose average length, height and thickness are respectively 28 cm, 16 cm and 8 cm. The liver is involved in many physiological functions such as the detoxification of metabolites, the secretion of hormones and proteins, the regulation of glycogen and the filtering of veinous blood to remove toxins.

Three types of blood vessels are located inside the liver: the hepatic artery bringing arterial blood (25% of the flow), the portal vein conveying blood from the digestive tract to be detoxified (75% of the flow), and the sushepatic veins collecting the filtered blood and draining into the inferior veina cava. The bile ducts carry the bile secreted inside the liver into the gallbladder through the cystic duct and to the pancreas through the common bile duct. The portal vein, hepatic artery and common bile duct jointly enter the liver through a central area called the hepatic hilum.

The liver is separated into a right and left lobes by the falciform ligament and the Couinaud classification further divides the liver into 8 segments. The liver consists of around one million lobules, hexagonal cylindrical structures of 0.8-2mm in diameter. Each lobule constitutes the basic functional unit of the liver and it is made of millions of hepatic cells (hepatocytes) connected to the arterial and the 2 veinous trees. The homogeneous spread of lobules and the amount of blood stored within the liver (around 450 mL of blood which amounts to 10% of the body's blood volume) explains that the mechanical behavior of the liver parenchyma can be assimilated at a certain scale to that of a porous medium. However, the global mechanical behavior of the liver is also impacted by the presence of its inner fibrous coat known as the Glisson capsule which ensheaths the veins, arteries and ducts within the organ.

In terms of mechanical boundary conditions, the liver is maintained in the abdomen through a two-layered fibrous membrane (peritoneum) that encases the liver in the abdomen except in a region where it connects directly to the diaphragm. Furthermore there are several ligaments among which the falciform and coronary ligaments connecting the liver to the diaphragm. The portal, arterial and biliary trees on one hand and the hepatic veinous tree on the other hand also play an important role in maintaining the liver in the abdominal cavity.

To conclude, the liver is a fairly complex organ to characterize mechanically due to its multiple components (lobules, blood, capsule) and its soft boundary conditions. The liver of a specific subject can be tridimensionnally reconstructed from CT scan or MR imaging. The injection of contrast agent during multiphase CT imaging allows to observe the arterial and venous trees and to detect the potential presence of lesions.

Finite Element models and constitutive laws of the liver: literature review

In addition to the definition of the geometric domain of the liver and its boundary conditions, it is necessary to define a suitable constitutive law in order to build a biomechanical finite element model of the liver. In this section, we review the most common non-linear constitutive laws to model the liver tissue behavior. This review is limited to the studies with model parameter identification from experimental tests. After referring to the main non-linear constitutive laws used for the liver tissue, the main physiological properties and their implications in the non-linear mechanical response are investigated.

Non-linear constitutive laws for liver modeling

The non-linear models dedicated to biological soft tissues are based on the identification of non-linear constitutive laws that are defined by a specific strain energy function W. In this subsection, the different forms of the strain energy functions identified from hepatic tissue tests are exposed and compared under uniaxial loading stress/strain relationship using

T = 2F •
∂W ∂C where T, F and C are the first Piola-Kirchhoff stress tensor (Lagrangian stress), the deformation gradient tensor and the right Cauchy-Green strain tensor, respectively. The most popular non-linear constitutive laws are summarized under the assumption of incompressibility in Table 1, (λ i , I i being the principal stretches and the invariants of the stretch tensor for i = 1, 2, 3 respectively).

Polynomial hyperelastic models

Polynomial forms of the strain energy are the most popular for the implementation of the liver tissue non-linear behavior. The Neo-Hookean model is the simplest polynomial function, based on the statistical thermodynamic and entropy conservation of cross-linked polymer chains. Hepatic tissue characterization using the Neo-Hookean model has been proposed in several studies [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF][START_REF] Ahn | Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations[END_REF]. The Mooney-Rivlin model consists of a generalization of the Neo-Hookean model, by the addition of the second invariant of the stretch tensor. This model is one of the most used at different orders to simulate the non-linear mechanical response of the hepatic tissue ( [START_REF] Hu | A biomechanical model of the liver for reality-based haptic feedback[END_REF] (N = 1), [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF] (N = 3), [START_REF] Umale | Experimental mechanical characterization of abdominal organs: liver, kidney & spleen[END_REF] (N = 2), [START_REF] Fu | Liver tissue characterization from uniaxial stress-strain data using probabilistic and inverse finite element methods[END_REF] (N = 2), [START_REF] Fu | Modelling and simulation of porcine liver tissue indentation using finite element method and uniaxial stress-strain data[END_REF] (N = 2)). However, the neo-Hookean model and the low order Mooney-Rivlin (up to the second order) alone have been shown to be not sufficient to reproduce the entire compression/elongation behavior of the hepatic tissue [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF]. The use of high order Mooney-Rivlin constitutive laws is thought to result in the identification of negative coefficients, which is not physically valid [START_REF] Hu | Characterization of soft-tissue material properties: large deformation analysis[END_REF]. Reduced polynomial models [START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF], Bogen [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF] and Ogden strain functions ( [START_REF] Hu | Characterization of soft-tissue material properties: large deformation analysis[END_REF] (N = 3), [START_REF] Gao | Estimating zero-strain states of very soft tissue under gravity loading using digital image correlation[END_REF] (N = 2 to 4), [START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF] (N = 4), [START_REF] Untaroiu | Material characterization of liver parenchyma using specimen-specific finite element models[END_REF] (N = 1), [START_REF] Lu | Effect of storage on tensile material properties of bovine liver[END_REF] (N = 1), [START_REF] Untaroiu | Modeling the biomechanical and injury response of human liver parenchyma under tensile loading[END_REF](N = 1)) have been also largely used, but there is still no clear consensus on which of the polynomial forms models the nonlinear behavior of the hepatic tissue the most accurately.

Exponential and logarithmic hyperelastic models

Exponential and logarithmic constitutive forms have been presented as potential models representing the stressstrain behavior of the hepatic tissue at least as well as the simple polynomial strain energy functions for both compression and elongation [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF]. Different formulations of exponential strain energy functions (such as the Fung-Demiray model) have been successfully identified on both the ex vivo and in vivo (mainly indentation tests) experimental response of the hepatic tissue [START_REF] Carter | Measurements and modelling of the compliance of human and porcine organs[END_REF][START_REF] Tamura | Mechanical characterization of porcine abdominal organs[END_REF][START_REF] Brown | In-vivo and in-situ compressive properties of porcine abdominal soft tissues[END_REF][START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF][START_REF] Roan | The nonlinear material properties of liver tissue determined from no-slip uniaxial compression experiments[END_REF][START_REF] Rosen | Biomechanical properties of abdominal organs in vivo and postmortem under compression loads[END_REF]. Both exponential and logarithmic models show similar capacities to simulate the hepatic tissue response [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF].

Combined models

Combining the different forms of the strain energy functions is thought to result in the improvement of the modeling of the hepatic tissue, including high-order polynomial formulations [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF]. The most popular approach consists in combining polynomial with either logarithmic or exponential constitutive laws [START_REF] Fu | Liver tissue characterization from uniaxial stress-strain data using probabilistic and inverse finite element methods[END_REF][START_REF] Fu | Modelling and simulation of porcine liver tissue indentation using finite element method and uniaxial stress-strain data[END_REF], with the possibility of including transverse isotropy using the fourth strain invariant [START_REF] Chui | Transversely isotropic properties of porcine liver tissue: experiments and constitutive modelling[END_REF]. The Veronda-Westmann strain energy function combines an exponential form with a polynomial formulation that includes the second strain invariant. Initially developed to simulate the skin response, this law has been successfully applied to healthy liver tissue [START_REF] Yin | Modeling of elastic modulus evolution of cirrhotic human liver[END_REF][START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF].

Inclusion of viscosity and porosity in the non-linear constitutive laws

By modification of the formulations, the influence of the viscosity has been included by some of the authors for both the Neo-Hookean model [START_REF] Miller | Constitutive modelling of abdominal organs[END_REF] and the first-order Mooney-Rivlin model [START_REF] Liu | Large deformations shear properties of liver tissue[END_REF][START_REF] Kim | Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse fe parameter estimation[END_REF][START_REF] Samur | A robotic indenter for minimally invasive measurement and characterization of soft tissue response[END_REF]. This extension of polynomial forms gives reasonable estimation of the hepatic tissue response at large strains over a large range of strain rates [START_REF] Liu | Large deformations shear properties of liver tissue[END_REF]. Rate dependence can also be added in a hyper-viscoelastic expression of the Ogden formulation by representing the relaxation function of the liver tissue as the Prony's series [START_REF] Sato | Hyper-viscoelastic response of perfused liver under dynamic compression and estimation of tissue strain thresholds with a liver finite element model[END_REF](N = 6). A non-linear strain-hardening fractional derivative constitutive law has been proposed for liver tissue by Nicolle et al. [START_REF] Nicolle | A strain-hardening bi-power law for the nonlinear behaviour of biological soft tissues[END_REF]. This approach gives the possibility to simulate the response of the tissue in both the linear and the non-linear viscoelastic regime over a significant range of compression strains (ranging from 0.01 to 1) and strain rates (ranging from 0.0151 s -1 to 0.7 s -1 ). Another form of constitutive laws consists of the combination of a polynomial, exponential or logarithmic non-linear strain energy function with viscosity by the incorporation of time-dependence in a simplified integral model [START_REF] Fung | Elasticity of soft tissues in simple elongation[END_REF]. While this Quasi-Linear Viscoelastic (QLV) formulation has been successfully used to model the liver [START_REF] Tamura | Mechanical characterization of porcine abdominal organs[END_REF][START_REF] Nava | Experimental observation and modelling of preconditioning in soft biological tissues[END_REF][START_REF] Roan | The nonlinear material properties of liver tissue determined from no-slip uniaxial compression experiments[END_REF][START_REF] Nava | In vivo mechanical characterization of human liver[END_REF][START_REF] Jordan | Constitutive modeling of porcine liver in indentation using 3d ultrasound imaging[END_REF], its capability to simulate the hepatic tissue at large strains and high strain rates has not yet been fully demonstrated.

In addition to viscosity, the inclusion of porosity in a second order reduced polynomial model has been proposed by [START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation[END_REF] where she uses biphasic theory to include flow-independent hyperelasticity [START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation[END_REF]. This approach consisting of considering the liver as a fluid-filled sponge has also been proposed by Raghunathan et al. in 2010 [40]. This model has then been extended by Marchesseau et al. in 2010 [41] and is further developed in this subsequent chapter. While playing a major role in the non-linear response of the liver, the inclusion of porous components in the constitutive hepatic models is still marginal.

Influence of the mechanical test conditions

Some experimental and physiological aspects play a major role in the capability of a constitutive law to model the mechanical response of the liver tissue under specific loadings. These characteristics have to be carefully defined and determine not only the choice of the constitutive law but also the choice of the experimental tests to identify the parameters of the rheological model.

From in vitro to in vivo characterizations

In most studies, the non-linear constitutive laws for the liver have been identified exclusively from ex vivo experiments. Owing to its invasiveness, ex vivo testing (cut samples or whole organ) is associated with possible limitations and biases (destructive invasive approaches) when compared to in vivo physiological states (pressurization, vascularization, moisturizing, post-mortem degradation and necrosis). The storage of the specimens is one of the first critical aspects of in vitro testing. It has been showed that the stiffness of the hepatic tissue increases not only after freezing/thawing (by a factor 1.44 to 1.68 at 20 % of strain after 30 days of frozen storage) but also with increased preservation time (by a factor 1.58 to 1.96 at 20 % of strain after 60 days of frozen storage) [START_REF] Lu | Effect of storage on tensile material properties of bovine liver[END_REF]. The local and global failure strains have been shown to decrease significantly with both freezing/thawing and increased storage time. More than 50 % and only 17 % differences have been observed in steady state stiffness by Ottensmeyer et al. between in vivo and unperfused liver tissue, and between in vivo and perfused hepatic tissue (that maintains temperature, hydration, and physiologic pressure), respectively [START_REF] Ottensmeyer | Tempest 1-d: an instrument for measuring solid organ soft tissue properties[END_REF][START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation[END_REF]. This observation illustrates the fact that perfusion is playing a major role in the non-linear response of the liver, which is often underestimated by the in vitro identifications. Moreover, the models identified by these different tests conditions are strongly modified in their time-domain and frequency domain responses, illustrating the importance of identifying the hepatic non-linear models on in vivo experimental measurements to account for viscosity in the liver hyperelastic models [START_REF] Jordan | Constitutive modeling of porcine liver in indentation using 3d ultrasound imaging[END_REF][START_REF] Sato | Hyper-viscoelastic response of perfused liver under dynamic compression and estimation of tissue strain thresholds with a liver finite element model[END_REF].

Up to now, the in vivo techniques used to identify the non-linear response of the liver tissue consisted either in indentation [START_REF] Carter | Measurements and modelling of the compliance of human and porcine organs[END_REF][START_REF] Ottensmeyer | Tempest 1-d: an instrument for measuring solid organ soft tissue properties[END_REF][START_REF] Hu | Characterization of soft-tissue material properties: large deformation analysis[END_REF][START_REF] Kim | Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse fe parameter estimation[END_REF][START_REF] Samur | A robotic indenter for minimally invasive measurement and characterization of soft tissue response[END_REF][START_REF] Jordan | Constitutive modeling of porcine liver in indentation using 3d ultrasound imaging[END_REF][START_REF] Ahn | Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations[END_REF][START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF] or aspiration [START_REF] Nava | Experimental observation and modelling of preconditioning in soft biological tissues[END_REF][START_REF] Nava | In vivo mechanical characterization of human liver[END_REF] or Monitorized Endoscopic Grasper [START_REF] Rosen | Biomechanical properties of abdominal organs in vivo and postmortem under compression loads[END_REF]. Although maintaining most of the in vivo parameters of the physiological environment, these techniques are limited to the local and superficial characterization of the organ and necessitate a direct contact with the liver, being so partially invasive. Magnetic resonance-and ultrasound imaging-based elasticity imaging (referred to as elastography) are the only methods giving the opportunity to identify in vivo and non-invasively the mechanical response of the organs [START_REF] Vappou | Magnetic resonance-and ultrasound imaging-based elasticity imaging methods: A review[END_REF]. The comparison between in vivo elastographic and ex vivo rheometric measurements has been investigated [START_REF] Chatelin | In vivo liver tissue mechanical properties by transient elastography: Comparison with dynamic mechanical analysis[END_REF]. While elastography is a continuously expanding field for diagnostic purpose, it is still limited to the investigation of the linear behavior of soft tissue and cannot yet been used for the non-invasive identification of the non-linear mechanical response of the liver tissue.

Dependence to the strain rate

Porosity and viscosity are two characteristics, which would both involve a time and frequency dependence in the response of the tissue. In most of the experiments aiming at identifying non-linear constitutive laws, these properties will be characterized by a strain rate dependence on the medium. By identifying non-linear models from uniaxial tensile/compression tests on ex vivo samples with varying loading rates (strain rates ranging from 0.003 to 0.606 s -1 ), the effect of strain rate on porcine liver has been shown to be relatively insignificant by Chui et al. [START_REF] Chui | Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation[END_REF]. However, by extending the range of strain rates to all the values proposed in the literature for the non-linear characterization of the hepatic tissue (ranging from 0.003 to 22.5 s -1 ), significant effects of the loading speed appear clearly, as reported in Figure 2(A) (non-linear models expressed in uniaxial tensile/compression as the first Piola Kirchhoff stress over the stretch ratio). These specificities have been confirmed from non-linear shear tests showing that the liver tissue has to be considered as a very soft solid or highly viscous fluid [START_REF] Liu | Large deformations shear properties of liver tissue[END_REF]. Tensile tests have shown that the elasticity increases from 58 to 61, 68 and 100 kPa with increasing strain rates from 0.01 to 0.1, 1 and 10 s -1 , respectively [START_REF] Lu | Effect of storage on tensile material properties of bovine liver[END_REF]. The strain rate dependence affects not only the non-linear behavior but also the failure stress with average 2nd Piola-Kirchhoff stress values ranging from 33 to 94 kPa for strain rates from 0.01 to 1 s -1 [START_REF] Untaroiu | Material characterization of liver parenchyma using specimen-specific finite element models[END_REF]. Due to variations in the protocols from the studies reported in Figure 2(A), it remains difficult to extract a quantitative description of the strain rate dependency from the literature. One of the easiest way to compare all of these results is to extract the equivalent elasticities at small strain through the chord modulus (equivalent linear behavior between 0 and 10 % strain) from this figure and to represent them as function of the strain rate for the different model categories (Figure 2(B) and (C) for the models identified from compression and tensile tests, respectively). The same trend is observed in both compression and tensile (named group 1 and group 2, respectively), with a similar stiffening with increased strain rate. Only one group (named group 3) shows lower strain rate dependence, but correspond to studies dealing with modeling of the rupture and nonlinear behavior at large strain only (with Ogdens models limited to the first order). The strain rate dependence shows the necessity of accounting for the strain rate influence in the non-linear liver models by means of the inclusion of porosity and viscosity (for the lowest and highest loading rates, respectively) in the constitutive law. The influence of the strain rate is fundamental and is conditioning the non-linear constitutive law used to simulate the liver mechanical response, depending on the application of the model (such as computer assisted surgery -low strain rates -or abdominal impact trauma -high strain rates) Figure 2: First Piola-Kirchhoff (nominal) stress as a function of stretch for the main studies dealing with the influence of the strain rate on the liver non-linear modeling in tensile/compression uniaxial testing (A). From this graph, the equivalent elasticity values are deduced as a function of the strain rate for the linear part (strains < 10%) in compression (B) and tensile (C) for the non-linear models proposed in the literature (color scale).

Anisotropy

Most of the biological soft tissues, such as muscle and brain, are composed of fibers, which results in highly anisotropic mechanical properties and behavior. Accounting anisotropy in the tissue model is a critical issue. However, divergences can be observed in the literature about the inclusion of anisotropy in the nonlinear models of the healthy hepatic tissue. While high transverse isotropy has been related in both tensile and compression by Chui et al. in 2007 (with a factor 2 between longitudinal and transverse stiffness at 20 % strain), the large majority of the studies observe and model the liver tissue as isotropic [START_REF] Chui | Transversely isotropic properties of porcine liver tissue: experiments and constitutive modelling[END_REF]. In 2011, no significant differences were found by Pervin et al. in the non-linear compressive response obtained from bovine specimens along and perpendicular to the liver surface for a large range of strain rates [START_REF] Pervin | Dynamic compressive response of bovine liver tissues[END_REF]. Similarly, most of the experimental results and models in the literature are indicating isotropic behavior of the liver tissue.

Inter-species variations

In spite of major morphological variations, most of the non-linear constitutive laws of hepatic tissue have been identified from animal specimens (rhesus monkey [START_REF] Miller | Constitutive modelling of abdominal organs[END_REF], swine [START_REF] Tamura | Mechanical characterization of porcine abdominal organs[END_REF][START_REF] Brown | In-vivo and in-situ compressive properties of porcine abdominal soft tissues[END_REF][START_REF] Sato | Hyper-viscoelastic response of perfused liver under dynamic compression and estimation of tissue strain thresholds with a liver finite element model[END_REF] and bovine [START_REF] Pervin | Dynamic compressive response of bovine liver tissues[END_REF]). Whatever the species is, the hepatic tissue has been shown to have highly non-linear stress-strain behavior and most of the current liver models include the use of animal tissue to simulate human organs. However, the disparity between the experimental protocols related in the literature does not allow to highlight significant interspecies differences in the non-linear behavior of the liver tissue. A comparison between the models identified using similar protocols by [START_REF] Untaroiu | Material characterization of liver parenchyma using specimen-specific finite element models[END_REF]2015 shows that the human liver tissue is 1.34 to 1.48 and 1.49 to 1.73 times stiffer than the porcine hepatic tissue at 10 % and 20 % strain, respectively [START_REF] Untaroiu | Material characterization of liver parenchyma using specimen-specific finite element models[END_REF][START_REF] Untaroiu | Modeling the biomechanical and injury response of human liver parenchyma under tensile loading[END_REF].

Characterization of the Glissons capsule

While no difference have been observed by testing hepatic tissue in presence and in absence of the capsule by some authors [START_REF] Chui | Transversely isotropic properties of porcine liver tissue: experiments and constitutive modelling[END_REF], the capsule has specific characteristics, which are thought to result in modification of the global nonlinear mechanical behavior of the liver. The addition of the capsule to the simulation process more accurately reflects the boundary conditions present in a real in vivo organ [START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF]. The capsule has been shown as being mechanically isotropic at a meso-scale. The capsule has a great influence on the fracture of the liver, involving ultimate strain between 47 and 50.5 % and ultimate load of 0.3 N/mm for the entire organ [START_REF] Brunon | Mechanical characterization of liver capsule through uniaxial quasi-static tensile tests until failure[END_REF][START_REF] Brunon | Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests[END_REF]. Based on identification of the third-order Ogden model and assuming incompressibility, isotropy and strain rate independence, the small strains and large strains elastic moduli have been obtained at 8.22 and 48.15 MPa by Umale et al., i.e. significantly stiffer than the hepatic tissue alone [START_REF] Umale | Experimental in vitro mechanical characterization of porcine glisson's capsule and hepatic veins[END_REF], confirming the importance of including the capsule to simulate the response of the entire liver.

Finite Element Models of the liver

The development of liver finite element models at the organ scale, and not only at tissue scale, has initially been driven by surgery simulation applications thus requiring near real-time computations. In this context, approaches often assumed a linear elastic behavior [START_REF] Cotin | A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation[END_REF] discretized on linear tetrahedron finite elements which naturally leads to solving a linear system of equations whose inverse could eventually be precomputed [START_REF] James | ArtDefo accurate real time deformable objects[END_REF].

Since linear elastic materials are not suitable for large displacements, several authors in the computer animation community have proposed corotational elastic models [START_REF] Suwelack | Quadratic Corotated Finite Elements for Real-Time Soft Tissue Registration[END_REF][START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF] where linear elastic stiffness matrices are rotated for each element. While coping with finite strains, those models are restricted by the material linearity. More realistic soft tissue deformations may be obtained by resorting to hyperelastic materials minimizing a continuum strain energy. For real-time computation, early approaches have been based on St Venant Kirchhoff materials [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF] which exhibit a linear stress-strain relationship. Significant speed-up can be obtained by using reduced basis of deformation [START_REF] Barbič | Real-time subspace integration for St. Venant-Kirchhoff deformable models[END_REF], Proper Generalized Decomposition [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a pgd approach[END_REF] or by grouping expressions on edges, triangles and tetrahedra [START_REF] Picinbono | Non-Linear Anisotropic Elasticity for Real-Time Surgery Simulation[END_REF] when discretized on linear tetrahedra.

For general hyperelastic materials, several authors have relied on the Finite Element Method to simulate soft-tissue deformation with explicit time integration schemes [START_REF] Speidel | Intraoperative surface reconstruction and biomechanical modeling for soft tissue registration[END_REF]. For instance, Miller et al. [START_REF] Miller | Total lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF] have developed the Total Lagrangian Explicit Dynamic (TLED) algorithm with Neo-Hookean materials. This approach has been combined with Prony series to model visco-elasticity and has been implemented on Graphics Processing Units (GPU) [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and gpu execution[END_REF] to reach real-time computations. However the main limitation of this approach is that it relies on explicit time integration schemes which greatly simplifies the update at each time step but requires small time steps to keep the computation stable especially for stiff materials.

A porous visco-hyperelastic Finite Element Model of the Liver

As previously discussed section 3, hepatic tissue material law is highly dependent of the strain rate for low and high deformations. The dependence in high deformation should be taken into account by adding a viscosity component to the hyperelasticity, while the dependence in low deformation is mainly due to porosity. We chose to combine a hyperelastic model of the hepatic tissue with viscosity and porosity as described in [START_REF] Marchesseau | Fast porous visco-hyperelastic soft tissue model for surgery simulation: Application to liver surgery[END_REF]. Figure 3(Left) shows schematically how these three components interact within a physically-based model. The model is made of an isochoric part containing the viscosity and hyperlasticity in series, and a volumetric part to account for the extracellular fluid present in the liver, acting in parallel to the isochoric component. The porosity component introduces the fluid pressure as an additional state variable. In the next sections, the individual components and underlying variables are described in details.

Optimised assembly of Finite Elements for Hyperelastic Materials

Limitations of the classical FEM approach

Finite Element models for hyperelasticity usually require significant computation times when solved by implicit time integration schemes. The objective of this section is to introduce a more optimised discretization method suitable for all hyperelastic materials that increases the computational efficiency when assembling force vectors and stiffness matrices, necessary to the resolution of the FEM. To discretize the liver geometry, we use tetrahedral linear finite elements built from 3D medical imaging. Linear tetrahedra have constant strain which implies using a single integration point and simplifies the computation of shape functions gradient. However, the developed optimizations could be easily extended to other elements such as linear hexahedra elements or high order tetrahedra.

As described in previous chapters, any hyperelastic material is fully determined by its strain energy function W which describes the amount of energy necessary to deform the material. Based on variational principles [START_REF] Zienkiewicz | The Finite Element Method[END_REF], classical approaches involve the expensive computation of the second Piola-Kirchhoff (SPK) stress tensor S and the elasticity tensor N i jkl = 2 ∂S i j ∂C kl where C is the right Cauchy-Green deformation tensor. This discretization method suffers three limitations:

• The formulation and the computation of the SPK stress tensor S and the elasticity tensor N can be fairly complex.

Indeed, the first and second derivatives of the Jacobian J = detF with respect to C is non trivial and involves the inversion of C. This makes the expression of the derivatives of the invariants of C, Ī1 and Ī2 particularly cumbersome and therefore computationally expensive to evaluate (see for instance [START_REF] Weiss | Finite element implementation of incompressible, transversely isotropic hyperelasticity[END_REF] for details).

• The strain-displacement matrix B involved in the force and stiffness calculations combines two terms: the deformation gradient F and the gradient of the shape functions D i . The former changes at each iteration while the latter is constant. For basic hyperelastic materials for which the elasticity tensor N is constant or nearly constant, this is not optimal and a better choice would be to isolate the deformation gradient and to combine the shape vectors with the elasticity tensor.

• The classical FEM formulation of hyperelastic material cannot cope with nearly flat (J ≡ 0) or even inverted tetrahedra (J < 0). Although such deformation are non physical and do not represent a meaningful configuration, in interactive simulations, such cases of extreme compression can be met due to contact with rigid objects or to non-physical user-defined gestures. For instance, if the user-interface is not equipped with force-feedback, the tool controlled by the user can undergo non plausible trajectories. To cope with this, it is important to propose an hyperelastic material which can handle nearly flat tetrahedra.

Multiplicative Jacobian Energy Decomposition (MJED)

Our original discretization method addresses at least partially the three limitations of classical FEM. It is important to note that the approach developed in this section is completely equivalent to the classical FEM one but leads to more efficient assembly of the matrices before resolution of the system. First, instead of computing the force ( f ) and stiffness matrix (K) using the first and second derivative of the energy with respect to C (leading respectively to S and N), we compute them directly by deriving the energy with respect to the nodal position (Q i ):

f i = - ∂W ∂Q i T and K i j = ∂ 2 W ∂Q j ∂Q i
Then, our approach consists of three separate contributions: i) Decomposition of strain energy We propose to write the strain energy functions as a sum of terms

W k = g k ( Ĩ)h k (J)
or a function of it (for instance its exponential for Fung's law), where Ĩ = (I 1 , I 2 , I 4 ...). The purpose is to decouple in the strain energy, the invariants of the Cauchy-deformation tensor C from the Jacobian J. This allows to avoid matrix inversions and complex derivative expressions. Therefore g k is independent of J, its derivative will not involve any matrix inversions. This decomposition applies to every studied cases (Costa's law, Veronda Westmann, Arruda-Boyce, StVenant Kirchhoff, NeoHookean, Ogden, Mooney Rivlin, and anisotropic versions of them). Using this decomposition of strain energy enables complex material formulation to be computed more efficiently with only a sum of reasonably simple terms and no matrix inversions. Once the decomposition is done, getting h k (J) = dh k dJ requires a 1D derivation, and getting

S k = 2 ∂g k ( Ĩ)
∂C requires to combine well-known derivatives of the invariants (such as

∂I 1 ∂C = Id or ∂I 2 ∂C = IdI 1 -C).
The nodal forces therefore only require the inputs of the gradient of the shape functions D i , and the strain energy terms h k , h k , g k and S k (see [START_REF] Marchesseau | Fast porous visco-hyperelastic soft tissue model for surgery simulation: Application to liver surgery[END_REF] for full derivation):

f i = -V 0 k       h k (J)g k ( Ĩ) ∂J ∂Q i T + h k (J)F S k D i      
ii) Formulation of the stiffness matrix Implicit time integration schemes require the computation of the tangent stiffness matrix at each time step. This naturally involves elasticity tensors computed as the derivative of S k for each tetrahedron and at each time step. MJED leads to far simpler expressions of those tensors because S k is independent of J. Furthermore, in many common materials, the term containing those elasticity tensors can be precomputed. The full expression of the stiffness matrix includes 6 terms:

∂ 2 W(T P ) ∂Q i ∂Q j = V 0 k (G k + H k + I k + Λ k + M k + R k ) where                  G k = ∂h k (J) ∂Q j T ∂J ∂Q i g k ( Ĩ) H k = h k (J) ∂ ∂Q j ∂J ∂Q i g k ( Ĩ) I k = h k (J) ∂g k ( Ĩ) ∂Q j T ∂J ∂Q i Λ k = ∂h k (J) ∂Q j T D T i S k F T M k = h k (J) D T i S k ∂F T ∂Q j T R k = h k (J) ∂S k ∂Q j D i T F T
We focus in more details on the term involving the third order tensor ∂S k ∂Q j :

R k = h k (J) ∂S k ∂Q j D i T F T
The third order tensor ∂S k ∂Q j is computed through the chain rule with the the elasticity tensor ∂S k ∂C . Instead of representing the fourth order tensor ∂S k ∂C through 6 × 6 matrices, we define it through its application ∂S k ∂C : H on any symmetric matrix H. More precisely, we define the scalars (β k l , δ k l ), and 3 × 3 symmetric matrices (A k l , B k l ) such that :

∂S k ∂C : H = l β k l A k l HA k l + δ k l (H : B k l )B k l
where A : B = tr(B T A) for any two matrices A, B. Finally, the term R k is found to be a combination of two terms:

h k (J)F L k l (i, j)F T and h k (J)F U k l (i, j)F T where L k l (i, j) and U k l (i, j) are linear matrices:

       L k l (i, j) = β k l A k l D i ⊗ D j A k l + A k l (D j • A k l D i ) U k l (i, j) = δ k l B k l D j ⊗ D i B k l
This formulation leads to an optimization for the assembly of the stiffness matrix for two reasons. First, only scalars and 3 × 3 symmetric matrices are involved in the computation. Second, except for the Ogden model, the matrices A k l and B k l are constant and therefore matrices L k l (i, j) and U k l (i, j) can be precomputed for each tetrahedron before the simulation.

iii) Coping with highly compressed elements In case of high compression, the volumetric terms h k (J) in the strain energy become dominant. This makes the stiffness matrix singular and thus leads to numerically unstable computations because there are an infinity number of deformed configurations leading to the same value of J. In order to cope with this, [START_REF] Teran | Creating and simulating skeletal muscle from the visible human data set[END_REF] perform an SVD decomposition of the deformation gradient matrix. To avoid this computationally expensive decomposition, we propose instead to regularize the term G k = h k (J) g k ( Ĩ) ∂J ∂Q j ⊗ ∂J ∂Q i by replacing it with the following expression :

G k = h k (J)g k ( Ĩ) (1 -) ∂J ∂Q j ⊗ ∂J ∂Q i + ∂J ∂Q j • ∂J ∂Q i Id
The closer is to 1 the closer the G k matrix is to a diagonal matrix. In practice, we set

= (1 -J) if 0 ≤ J ≤ 1, = 0 if J ≥ 1 and = 1 if J ≤ 0.
In all cases, the trace of the regularized matrix is equal to the trace of the original matrix. By only regularizing the stiffness matrix, we still minimize the strain energy and therefore do not change the nature of the hyperelastic material. With this technique, it is even possible to handle inverted elements when the strain energy remains finite as J = 0. One limitation of the MJED method is that is does not directly compute the (first or) second Piola Kirchhoff stress tensor S = 2 W ∂C . If such stress tensors are required, then they need to be computed at extra cost (see below).

Visco-hyperelasticity based on Prony series

To model accurately the viscoelasticity of the liver, we propose to rely on Prony series [START_REF] Taylor | On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and gpu execution[END_REF]. This consists in adding to the hyperelastic SPK stress tensor now denoted h S some time dependent stresses. This time dependence is given by α(t) = α ∞ + i α i e -t/τ i with the condition α ∞ + i α i = 1. The visco-hyperelastic SPK tensor noted v S can then be written as:

v S = t 0 α(t -t ) ∂ h S ∂t dt = h S -i γ i where γ i = t 0 α i (1 -e (t -t)/τ i ) ∂ h S
∂t dt After a discretization over time this results in the recursive formula between time n -1 and time n:

γ n i = a i h S n + b i γ n-1 i where a i = ∆tα i ∆t + τ i and b i = τ i ∆t + τ i
∆t is the time step used for discretization and has to be the same as the time step for any solvers during the simulation. To combine the Prony series with our optimized hyperelastic formulation, we therefore need to compute the total hyperelastic SPK stress tensor h S n . This is done for each time step n (dropping the index for clarity) by computing the inverse deformation gradient :

h S = F -1 k ( f k (J)g k ( Ĩ)JF -T + f k (J)F h S k ) where F -1 = 4 l=1 P l ⊗ ∂J ∂Q l /J
The visco-hyperelastic nodal forces v f are therefore related to the hyperelastic ones h f by

v f i = h f i + V 0 F i γ i D i        D i Id
Adding the viscous properties through the Prony series does not have a significant impact on the total computation times despite the evaluation of the time dependent stresses γ n i and F -1 .

Poro-elasticity

We follow Kerdok's porosity model [START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation[END_REF] and consider the liver as a fluid-filled sponge. This model is made of a volumetric component (represented by σ Heq Figure 3) and a fluid phase with variable P f luid ). The volumetric component is governed by Hencky's elasticity [START_REF] Xiao | Hencky's elasticity model and linear stress-strain relations in isotropic finite hyperelasticity[END_REF] in which the Cauchy stress depends on the proportion of free-fluid (e.g. blood, water) in the liver parenchyma in the reference configuration (r w ) and the effective volumetric Jacobian J * = (r w + J -1)/r w by: σ Heq = K 0 r w ln(J * ) where K 0 is the bulk modulus of the material. With this model, when J get close to 1r w , the solid phase of the liver is completely compressed and the resulting stress is infinite. To avoid instabilities due to this infinite stress, we substitute σ Heq when J ≤ J 0 by its tangent curve at J 0 (see Figure 3(Right)). We set J 0 = 1r w + K 0 /K lim where K lim is a bulk modulus and represents the slope of the tangent. The fluid phase of the liver also applies some volumetric stresses due to the transient response of the fluid through the porous liver parenchyma. A straightforward way of modeling the porous behavior is through the linear Darcy's law. In this setting, variation of fluid pressure P f luid is governed by the variation of volume change and a diffusive process:

1

K lim Ṗ f luid = κ∇ 2 P f luid - J J
where κ is the permeability parameter. In Kerdok's model, the permeability κ is a function of J, but we propose to keep it constant to decrease its computational cost. Finally, the total Cauchy stress response in the volumetric part is defined by summing the solid and the fluid terms:

σ p = σ heq Id -P f luid Id
The Cauchy stress is translated as a poro-elastic force:

p f i = -σ p ∂J ∂Q i T V 0
added to the visco-hyperelastic forces. The additional stiffness matrix due to the porous behavior, which is regularized similarly to previously described in section 4.1.2 iii), is given as:

p K i j = V 0       ∂σ p ∂J ∂J ∂Q j T ∂J ∂Q i + σ p ∂ 2 J ∂Q j ∂Q i      

Results

MJED validation

First, we compared the computation time between the MJED and the classical FEM method, referred to as "Standard FEM", implemented in SOFA 1 . The results are given in Figure 4. We measured the time elapsed for the computation of the nodal forces and the stiffness matrices averaged over 100 iterations. For a more detailed differentiation between force and stiffness assembly computation, the reader is referred to [START_REF] Marchesseau | Fast porous visco-hyperelastic soft tissue model for surgery simulation: Application to liver surgery[END_REF]. We simulated the deformation of a cube with 20 700 tetrahedra and 4300 nodes. For all modes implemented the proposed strategy is definitely more efficient than the standard FEM, up to five times as fast for St Venant Kirchhoff material.

Second, the accuracy of the MJED computation was compared against the Open Source software FEBio (version 1.1.7) 2 where several elastic and hyperelastic materials are implemented and against analytical solutions for simple cube model. The mean relative difference is around 10 -6 for every models tested. Moreover, a mesh convergence study was performed to evaluate the trade-off between accuracy and speed on several meshes. Dividing a mesh 20 times only increases the accuracy by 1% while multiplying the computation time by 10.

Rheological testing on porcine liver: experimental materials and methods

To calibrate the visco-elastic parameters of our liver model, tests were performed on porcine livers. Dynamic viscoelastic behavior of hepatic tissue was investigated using in vitro Dynamic Mechanical Analysis (DMA) in rotating shear. Tests were carried out on liver cylindrical shaped samples coming from five adult pigs (weighting between 25 and 35kg). Immediately after hepatectomy performed on anesthetized animals, entire livers were stored in an insulated container at 6C surrounded by ice. Cylindrical hepatic samples of 20 mm diameter and 41 mm thick were cut and tested within 6 hours post-mortem time. To avoid mechanical difference due to samples localization, 4 samples were tested for each of right, middle and left liver lobe. Hepatic tissue is considered as isotropic. At least 60 samples were tested (5 animals x 3 lobes x 4 samples).

Dynamic Strain Sweep tests and Dynamic Frequency Sweep tests were performed on a dedicated stress-controlled AR2000 (TA-Instruments, New Castle, DE, USA) rheometer in a parallel-plate configuration. A pre-compression of 5mN was applied and sand paper was fixed to the rheometer plates to insure grip with tested sample.

-Dynamic Strain Sweep

Those tests aim at giving linearity limit of the material elasticity. The sweep covers strain range from 0.01% to 20% which is sufficient according to the literature. Measurement of both storage G' and loss G" shear moduli as functions of shear strain showed that the linearity domain extends up to 1% to 2% of deformation.

-Dynamic Frequency Sweep These experiments were carried out in the linear viscoelastic strain range of the sample (gamma0=0.1%) with frequencies ranging from 0.1 to 4 Hz. From the 60 tested samples, we display the mean curves (Figure 5, mean +/standard deviation values) of the storage G' and loss G" shear moduli, which correspond to the real and imaginary parts of the complex shear modulus G * (ω) = G (ω) + iG (ω), respectively. From these results the Dynamic Modulus G can be obtained as a function of the frequency or function of the time, and the viscoelastic behavior can be modeled after fitting a generalized Maxwell model with two modes of relaxation to those measurements:

G(t) = G 0 (g ∞ + g 1 e -t/τ 1 + g 2 e -t/τ 2 )
where G ∞ = G 0 g ∞ is the equilibrium modulus, g 1 , g 2 , τ 1 , τ 2 are model relaxation parameters such as g ∞ + g 1 + g 2 = 1. The parameters are given Table 2. Our estimated initial and equilibrium shear moduli were 770 Pa and 333 Pa, respectively. While corresponding to one of the lowest values related in the literature, these parameters are in agreement with most of the studies dealing with hepatic tissue rheological testing at surgical strain rates. [START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF][START_REF] Tamura | Mechanical characterization of porcine abdominal organs[END_REF][START_REF] Kim | Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse fe parameter estimation[END_REF][START_REF] Ahn | Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations[END_REF][START_REF] Fu | Modelling and simulation of porcine liver tissue indentation using finite element method and uniaxial stress-strain data[END_REF] From the rheological experiments described in the previous section, we derive the initial shear modulus G 0 required in the hyperelastic term (e.g. Arruda-Boyce material) and the Prony series parameters required in the viscous term. To check that those parameters are indeed valid, we propose in this section to compare in silico simulations with the performed in vitro rheological tests. First, we check that the linearity domain for our hyperelastic materials matches the ones observed in the dynamic strain sweep experiments. We simulate a cylinder in extension in SOFA for several longitudinal stresses and estimate the associated strains. The linearity limit corresponds to the one given by experiments (1% -2%).

Second, dynamic frequency sweep tests have been simulated using similar geometries and boundary conditions than the DMA tests. An oscillating torque (amplitude M) is applied on a small cylinder (radius r = 10mm, height h = 4mm) at various frequencies ω. The amplitude of the torque is chosen so as to stay in the linear domain. The angle of rotation θ of the cylinder is measured as a function of time. This angle describes a sinusoidal curve which follows the torque amplitude with a shifted phase δ. Specific constraint is applied on the top cylinder nodes to enforce a pure rotation of those nodes (as to reproduce the pure grip of sand paper).

Using similar calculation of stress and strain as for the rheometry experiments, we can estimate the values of the Storage and Loss moduli to be compared with experimental data. It can be seen in Figure 5 that the simulation manages to capture the viscous behavior of the liver for small deformations with a mean relative error of 5%. By comparison with the constitutive laws in both compression and tension from the literature, this behavior observed in shear is in good agreement with most of the studies at similar strain rates (ranging from 0.008 to 0.250 s-1) [START_REF] Lister | Development of in vivo constitutive models for liver: application to surgical simulation[END_REF][START_REF] Tamura | Mechanical characterization of porcine abdominal organs[END_REF][START_REF] Kim | Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse fe parameter estimation[END_REF][START_REF] Ahn | Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations[END_REF][START_REF] Fu | Modelling and simulation of porcine liver tissue indentation using finite element method and uniaxial stress-strain data[END_REF]. By the use of strain rates similar to those from most of the surgical applications, this results confirms the capability of our model to simulate the liver response during surgical interventions.

20

Table 2: Values of the parameters identified from the DMA Tests for the viscosity and from the literature for the porosity.

Liver Simulations

To describe the influence of each component in the complete model several simulations were performed where the hyperelastic material used was a fifth-order Arruda-Boyce. The liver mesh was segmented from a CT image and meshed with the GHS3D software (1240 vertices and 5000 tetrahedra). An Euler implicit time integration scheme is used and the linear equations are solved with a conjugated gradient algorithm. As boundary conditions, several nodes of the liver are fixed along the vena cava and suspensive ligament. The liver deforms under the action of gravity such as to overpass the linearity limit of the material. All computations were performed on a laptop PC with a Intel Core Duo processor at 2.80 GHz. (i) Influence of the viscous component Adding viscosity to hyperelasticity increases the amplitude of the oscillations as the material becomes less stiff. Contrary to essentially hyperelastic model, the final state is really different from the initial state (see Figure 6 (Bottom)). Indeed, the use of Prony series leads to a multiplication of the SPK tensor by 1a k at infinte time. The frame rate is around 9 FPS against 10 FPS for hyperelasticity alone. However the implicit integration scheme allows larger time step (0.3s for instance) which makes the real-time interaction possible. High amount of extension and compression are possible which may be somewhat unrealistic, therefore the porous component is necessary to control the amount of viscosity. (ii) Poro-hyperelastic simulation We have implemented the porous component in parallel to the hyperelastic component using parameters based on Kerdok's [START_REF] Kerdok | Characterizing the nonlinear mechanical response of liver to surgical manipulation[END_REF] experimental data, shown Table 2. The simulated fluid pressure field during the deformation is shown in Figure 6 (Top) as a color map, ranging from dark blue (initial pressure) to red (highest pressure). Highest pressure in the fluid occurs when the liver is compressed either by gravity (diffusion starts at the top), either by elastic reaction (diffusion starts from the bottom).

(iii) Complete model Adding porosity to visco-hyperelasticity prevents the liver from having unrealistic large deformations. The deformation is no longer isotropic and changes over time. The addition of this component decreases the computational efficiency (6 FPS) since a semi-implicit integration scheme is used for the porous component. However, because of the fast variation of the explicit term J/J, the time step has to be decreased to 0.15s. On our laptop PC, the simulation is still fluid enough to allow user interactions.

Discussion and Perspectives

Liver biomechanics has been quite extensively studied in prior work and this chapter provides a synthetic view of hyperelastic material parameters that were fit on various rheological experiments. It appears that strain rate has a major influence on the behavior of the tissue behavior thus hinting for the use of visco-hyperelastic or even viscoporo-hyperelastic models. However, the variety of experimental conditions (ex-vivo vs in-vivo, strain rate, considered species . . .) makes it very difficult to provide a precise quantitative characterization of the liver mechanics. Clearly, significant advances in this direction can only be obtained by a coordinated effort by the scientific community to standardize the experimental process and the modeling choices.

Furthermore, the liver mechanical properties are strongly influenced by the physiological conditions of the subject. For instance, liver fibrosis is fairly common in chronic liver diseases and the hepatic stiffness has been shown to increase with the staging liver fibrosis and with the METAVIR scoring system [START_REF] Castra | Prospective comparison of transient elastography, fibrotest, apri, and liver biopsy for the assessment of fibrosis in chronic hepatitis c[END_REF] due to the excessive accumulation of collagen fibers. Fatty liver disease (hepatic steatosis) also translates into a significant increase in liver stiffness [START_REF] Huwart | Magnetic resonance elastography for the noninvasive staging of liver fibrosis[END_REF]. Elastography imaging techniques based on Ultrasound [START_REF] Gennisson | Ultrasound elastography: principles and techniques[END_REF][START_REF] Fatemi | Ultrasound-stimulated vibro-acoustic spectrography[END_REF] or Magnetic Resonance Imaging [START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propagating acoustic strain waves[END_REF] are currently being developed to provide patient specific maps of linear mechanical parameters (mainly elastic or viscoelastic shear moduli) of the main abdominal organs. Those non-invasive measurements are currently limited to small deformations but open new avenues for the design of patient specific biomechanical models. In particular, they could be useful to provide quantitative information about the boundary conditions of the liver with respect to neighboring organs for which little is known.

Finally, global finite element models of the liver have been mainly developed for medical applications and in particular for surgery training and therapy guidance. In both cases, fast or real time simulations are required thus leading to fairly simplified mechanical models. In this paper, a porous visco-hyperelastic model was presented with a fast assembly of the stiffness matrix leading to optimized performances. This approach can be further improved by taking into account the mechanical resistance due to the large vessels inside the liver and more appropriate boundary conditions. In addition, with the development of non invasive elastography and MR or ultrasound intraoperative imaging, it is foreseeable that such mechanical models could be validated and personalized with regional parameters and not only global ones.

Figure 1 :

 1 Figure 1: Liver Anatomy including the right and left lobes separated by the falciform ligament

Figure 3 :

 3 Figure 3: (Left) Representation of the constitutive model combining viscosity, hyperelasticity and porosity. (Right) Representation of the static Cauchy stress before and after substitution. (Red) Initial stress, (Green) Tangent curve, (Blue dots) New stress. Here f w = 0.8.

Figure 4 :

 4 Figure 4: Comparison of the computation times of nodal forces and stiffness matrices between two different discretization methods averaged over 100 iterations.

5. 3 .Figure 5 :

 35 Figure 5: Comparison of the simulated values with the data obtained by DMA testing. The moduli are given on a log-log scale. The material is St Venant Kirchhoff, similar values are found for other materials.

  Pa) g 1 (Pa) τ 1 (s) g 2 (Pa) τ 2 (s) r w K 0 (Pa) K lim (kPa) κ(m 4 /N s)

Figure 6 :

 6 Figure 6: (Top) Pressure field of the porous component on a liver under gravity, during one oscillation. (Bottom) Addition of viscosity to hyperelasticity: Comparison of the maximum amplitudes and final states. (Black) Initial position, (Blue) Hyperelastic liver, (Pink) Visco-hyperelastic liver.

Table

SOFA is an Open Source medical simulation software available at www.sofa-framework.org

FEBio is an opensource software package for FE analysis available at mrl.sci.utah.edu


Renaud Larrieu

Catherine Quilliet

Aurélie Dupont

Philippe Peyla

Collective orientation of an immobile fish school, effect on rheotaxis

We study the orientational order of an immobile fish school. Starting from the second Newton's law, we show that the inertial dynamics of orientations is ruled by an Ornstein-Uhlenbeck process. This process describes the dynamics of alignment between neighboring fish in a shoal -a dynamics already used in the literature for mobile fish schools. Firstly, in a fluid at rest, we calculate the global polarization (i.e. the mean orientation of the fish) which decreases rapidly as a function of the noise. We show that the faster a fish is able to reorient itself, the more the school can afford to reorder itself for important noise values. Secondly, in the presence of a stream, each fish tends to orient itself and swims against the flow: the so-called rheotaxis. So even in the presence of a flow, it results in an immobile fish school. By adding an individual rheotaxis effect to alignment interaction between fish, we show that in a noisy environment, individual rheotaxis is enhanced by alignment interactions between fish.

INTRODUCTION

The appearance of self-organization within a group of active entities is a fascinating phenomenon [1]. It has been studied for micro-organisms [2,3] as well as for active synthetic particles [4,5] and at larger scales for animals [6]. Fascinating self-organisation is also observed within fish shoals [START_REF] Radakov | Schooling in the Ecology of Fish[END_REF][START_REF] Pavlov | [END_REF][9][10]. Milling and schooling are collective phenomena occuring on scales much larger than an individual fish. This phenomenon has been studied since the begining of the 20 th century [11]. The structure and the sensitivity to external factors such as water temperature, light and darkness were analyzed by Breder [12]. Individual fish in a school were observed to swim for a longer duration when aligned, with lower tail-beat frequencies, smaller energy dissipation and respiratory rates, compared to fish swimming alone [13][START_REF] Ross | American shad[END_REF][15]. In addition, shoaling and alignment between fish are established as a result of many social and sensory factors like metabolism [16], alignment by vision [17] or food [18]. Recently, the study of out-ofequilibrium active systems [19] allowed scientists to substantially improve their knowledge in modeling this remarkable phenomenon. In the seminal work of Vicsek [20], an individual animal (bird or fish) adopts instantaneously the average orientation of its neighbors in the group, resulting in a collective motion that can be destroyed by noise. The noise source can be intrinsic to the fish or due to external conditions such as turbulent fluid flow [9]. Since then, more sophisticated force models have emerged that reproduce quite well the real behavior of schools of fish [21][22][23]. That class of social model allows to study several situations with some flexibility [24].

If collective motions have been extensively studied [1], quite poor literature is devoted to immobile groups of fish [25] which stay at the same place relatively to their living environment. Immobile fish schools can be observed in various situations and especially in reef regions [START_REF] Keenleyside | Diversity and Adaptation in Fish Behaviour[END_REF]. Origins of such an immobile state are diverse. It is likely that schools of fish that stop their movements and remain motionless for a period of time may achieve perceptual benefits [10]. Simultaneous stopping of a school of fish provides relatively quiet intervals to allow reception of potentially critical environmental signals, fish under predator threat that form non-moving "look around shoals" [START_REF] Radakov | Schooling in the Ecology of Fish[END_REF] may be an example. However, the most frequent origin of immobile school is a rheotactic effect that allows the fish to orient against a stream [START_REF] Pavlov | [END_REF] and is the object of the present model. This effect was studied in details by Potts [START_REF] Potts | [END_REF]. A school of the snapper Lutjanus monostigma was observed during several days and self-organized into a polarized and immobile school when submitted to tidal flow. Each fish were heading into the current in order to maintain their position by positive rheotaxis. This is done by swimming gently at an equal and opposite speed to the current. Indeed, by pointing ahead in a direction opposed to the flow can help the school to maintain its immobility in a region where food is present. A fish can individually find the direction of flow through sentitive captors [28,29] and can also try to align with its congeners.

In the following, we will first present the model of fish orientation with respect to neighbours and flow. We then show that alignment interactions within a shoal can increase rheotaxis efficency of a single fish.

MODEL

Let's consider a motionless fish (fig. 1) located at a fixed position and with a time-varying orientation, living in a school and interacting with its neighbors while attempting to orient itself in a direction opposite to a uniform flow (rheotaxis). Here, for simplicity we use a 2D fish shoal [22,23] with a circular shape of radius R. Each of the N fish is a discoid of radius a. When we vary N , we maintain the density ρ = N (a/R) 2 constant, typically, ρ ≈ 0.5. This corresponds to a quite dense school which can be oftenly encontered [START_REF] Simmonds | Fisheries Acoustics: Theory and Practices[END_REF] and which allows us to have small fluctuations of the local density of fish. In the spirit of social force model originally developped by Helbing [START_REF] Helbing | [END_REF]32], we consider here social torques. These "torques" are a measure for the internal motivations of each individual to perform certain movements (rotation) depending on its environment. Each fish needs to orient itself in the same direction as its neighbors within a chosen radius of 5a (i.e. 2.5 times its own size) and against the flow. The typical size 5a has been chosen to capture neighbors that are in the close neighboring of a given fish regarding the chosen density ρ. Starting with the second Newton's law for rotating bodies we can write:

I ωi = -ζω i + j∈Vi T I ij (φ ij ) + T R i (φ i ) + η(r i , t), (1)
FIG. 1. Fish in a flow. The flow is oriented toward the y < 0-direction. Angles: φi is the angle between orientation ei of the i th fish and the horizontal x-axis; θij is the angle between fish i orientation ei and the fish i-j orientation vector eij; finally φij is the angle between orientaions of fish i and j.

where I is the moment of inertia of any fish (supposed identical), ω i is the angular velocity of fish i located in r i and oriented with an angle φ i with the x-axis at time t. The angular acceleration is ωi .

Fluid friction is ζ. The i th fish interacts with its close congeners indexed by j in the V i ensemble of its N i neighbors (i.e. within a circle of radius 5a around fish i). This interaction is represented by a social torque T I ij which depends on the relative orientation between fish i and j: φ ij = φ i -φ j . The torque associated with the rheotaxis is T R i which depends on the orientation of the fish i with the y-axis chosen as the direction of the incident and uniform flow. Finally, the dynamics is perturbed by a noise term η(r, t) with < η(r, t) >= 0 and < η(r, t)η(r , t ) >∼ δ(r -r )δ(t -t ). Note that in our model, we do not consider an interaction term that depends explicitly on fish interdistance. This equation can be rewritten as:

ωi = - 1 τ ω i + 1 τ ω * i (φ ij ) + η(r i , t)/I, (2)
which expresses that each fish adjusts its angular velocity ω(t) towards a time-dependent target value ω * i = ω I * i + ω R * i , depending both on the fish-fish interaction (ω I * i ) and on the rheotaxis (ω R * i ) within an external flow. Expressions of ω I * i and ω R * i are given below. We use the time τ associated with dissipation as the characteristic time and we rescale by τ the other times associated with alignments and rheotaxis. Note that if we consider a fish with a size around a few cm, the time τ associated with dissipation during a solid rotation is about a few seconds which is much bigger than a typical time of reaction for alignment closer to a few tenths of a second. However, the time associated with dissipation can be much shorter if we consider that usually a fish is flexible and a change of orientation is driven by a deformation of its body which can reduces very much its inertia [33]. Following the spirit of [22,23] in order to describe the interaction between close fish, we write :

ω I * i = I τ 1 N i j∈Vi sin( φ ij 2 ) 1 + cos(θ ij ) 2 , ( 3
)
where I is the dimensionless amplitude of alignment interaction. The term sin(φ ij /2) accounts for alignment between fish i and j. If fish i and j are aligned in the same direction, it reads sin(φ ij /2) = 0 and then ω * i = 0 since fish i has no reason to rotate. But if fish i and j are anti-parallel, i.e. φ ij = ±π, then sin(φ ij /2) = ±1, since fish i must rotate. The term (1 + cos θ ij ) /2 is designed to ensure a frontal preference and some kind of rear blind angle [23]. In order to model the alignment against the flow (rheotaxis), we have:

ω R * i = F τ sin π/2 + φ i 2 . (4)
The term F is dimensionless and represents the amplitude of the rheotaxis. If φ i = -π/2, the fish does not rotate since it is already aligned against the flow (i.e. pointing in the y > 0 direction) it reads ω R * i = 0. But if the fish is aligned along the flow, (φ i = +π/2) the fish must turn back in order to point against the flow with the target angular velocity ω R * i = F/τ . Using an Euler-Maruyama integration [34], eq.( 2) reads:

dω i (t) = - dt τ [ω i (t) -ω * i (t)] + σN √ dt, (5)
also known as an Ornstein-Uhlenbeck process [21,22]. The noise amplitude is σ and N is a random gaussian variable of mean 0 and variance 1. For large values of time, this equation becomes stationnary (see appendix). In the absence of rheotaxy (F = 0), we can rescale time and angular velocities by τ /I and I/τ respectively. At stationnarity, it is easy to show that there is only one dimensionless number σ = στ 3/2 /I 1/2 that compares the amplitude of noise and the amplitude of alignment interaction (see appendix).

In the presence of rheotaxis (F = 0), three terms should be compared: the alignment interaction (amplitude I), the rheotaxis (amplitude F ) and the noise term (amplitude σ). Since we would like to vary F at constant I, we choose to rescale time and angular velocities by τ /F and F/τ respectively. Then, we get two dimensionless numbers I/F and F 1/2 /(στ 3/2 ). The last term being equal to (F/I) 1/2 /σ.

In the following, we will integrate numerically eq.( 5) in the absence or in presence of rheotaxis. For each situation, we will plot the global polarization (defined below) as a function of the above dimensionless numbers.

RESULTS.

We first analyse the alignment within an immobile school of N = 200 fish as a function of the noise σ and without rheotaxis (F = 0). Note that we tried different numbers of fish: N = 400, N = 200, We compute the mean value of the global polarization [20,22] defined as P = 1 N N i=1 e i :

P = 1 N N i=1 cos φ i (t) 2 + N i=1 sin φ i (t) 2 . (6)
The global polarization of the school is P = 1 when all fish point in the same direction while P = 0 means that the fish point in different directions.

The school initial orientation is polarized in a given random direction. After a transient time of order τ /I, fish can loose partially their mutual alignment because of the noise. Thus, we measure P for large values of time (t τ /I) as a function of σ (see fig. 2). We found that P drops abrupltly for σ ∼ 0.27. Here, the rescaled noise being σ = στ 3/2 / √ I, it means that, fast reacting fish (small values of τ /I) are more able to line up in a stronger noisy environment (large value of σ) than fish with a larger τ /I value that cannot sustain the same amplitude of noise to form a polarized school. Note also that when P = 0.6 ± 1, we obtain a maximum of polarization fluctuations whatever N (see appendix -fig. 8).

To visualize the loss of orientation when noise increases, let us first consider the map of individual polarizations around σ = 0.27. In fig. 3, we plot the time averaged value of each direction ēi (t) = 1/t max tmax 0 e i dt for different values of σ. We use t max /(τ /I) = 10 4 and dt/(τ /I) = 10 -2 . Starting with σ = 0.18, the group of fish is well polarized (fig. 3.a). Then, we increase the rescaled noise value. We see that for σ = 0.27, (fig. 3.b) a weak polarized region appears. Around σ = 0.54 (fig. 3.c) several weak polarized zones have invaded the school, leading to islands of polarized fish separated by unpolarized zones. For larger values σ = 1.0 (fig. 3.d), an entire unpolarized fish school (a so-called shoal) remains. We also calculate the correlation function C = e i .e j = 1/t max tmax 0 e i .e j dt for different distances d between fish i and j. As shown in figure 4, C decreases exponentially as a function of the fish-to-fish distance with a typical correlation length λ. This correlation length decreases by increasing the noise (fig. 5). We observe a small plateau for λ close to the school size R probably due to a boundary effect.

To study the effect of the alignment interaction between fish on the rheotaxis of the whole group, we now consider non-zero values of F . In the absence of alignment interactions between individuals (I = 0), each fish tends to orient itself against the flow (pointing toward y > 0 direction). The presence of noise perturbs the rheotactic behavior of each fish and the polarization drops by increasing the noise σ or τ and decreasing the rheotaxis F (fig. 6). We assume that even if a given fish is not perfectly oriented against the flow it still maintains its position within the shoal in order to stay with its congeners. Now, by switching on the alignment interaction between fish (i.e. for non-zero and positive values of I), we observe a clear increase of the global polarization (see figure 6 (black arrow)). Note that the x-axis ( √ F /(στ 3/2 )) is inversely proportional to σ. In the inset of fig. 6, we plot the polarization difference ∆P = P (I/F ) -P (I/F = 0) between the global polarization P of the school in presence of interactions and in the absence of social interaction (I = 0). We see that for large values of I/F , a strong increase of the polarization against the flow is obtained and reaches a maximum around F 1/2 /(στ 3/2 ) ≈ 1. For small rheotaxis or strong noise the global polarization drops to zero since the fish are pointing in all directions. On the contrary, for large rheotaxis or small noise all the fish are pointing in the direction opposed to the flow and P saturates. In both cases, the role of fish-fish interaction is inefficient. But between these two extreme cases (F 1/2 /(στ 3/2 ) ≈ 1), we observe a maximum of ∆P corresponding to a significant gain of rheotaxy by the interplay of the fish to fish interactions. It can be concluded that a strong fish to fish interaction improves the collective rheotaxis even if the individual rheotaxis is weak.

CONCLUSION

In this paper, we have studied the collective orientation of an immobile group of fish with two ingredients: a social torque to align fish with their close neighbors and an environmental torque to align fish with an external flow. We have modeled the inertial dynamics of groups of fish in the presence of noise. In the absence of an external flow, we show that for large values of a dimensionless noise σ 1, the group cannot globally polarize. In the presence of a flow, we show that strong social interactions help the group to detect and align even with weak individual rheotaxis. This model can be extended to the case of moving fish which for certain species are able to detect low gradients of velocities (in a Poiseuille flow) [35]. In this case, it would be interesting to study whether collective social interactions can improve the efficiency of this specific rheotactic behavior. Despite the simplicity of the model which does not account for hydrodynamic drag, we believe that several dynamics of group organization (or collective behaviors) in a complex environment can be captured.

where t = tF/τ and ω = ωτ /F . The target is:

ω * i = I F 1 N i j∈V sin( φ ij 2 ) 1 + cos(θ ij ) 2 -sin( π/2 + φ i 2 ), (10)
In this case, two dimensionless numbers are to be considered I/F and √ F /(στ 3/2 ).

Stationnarity

For large values of time ( t 1), we can assume stationnarity. In the presence of rheotaxis, if this hypothesis is true, we should have:

[ω i -ω * i ] dt = στ 3/2 F 1/2 N dt. ( 11
)
But this equality is not easy to prove numerically because of the presence of noise. So let us average each member of eq. ( 11) on time and fish. By integrating on time and using Ito isometry we obtain the following equality:

1 tmax tmax 0 [ω -ω * ] dt 2 = στ 3/2 F 1/2 2 ( 12
)
FIG. 9. Term T1 (lhs of eq.( 12)) as a function of T2 (rhs of eq.( 12)), Data are circles and the red solid line represents the first bisector, showing that T1 = T2.

with tmax = t max /(τ /F ) = 10 . In figure 9, we plot T 1 as a function of T 2 for different values of σ, τ , I and F . We show that these two terms are identical and thus stationarity hypothesis is true when tmax 1).

FIG. 2 .

 2 FIG. 2. Decrease of the global polarization as a function of the dimensionless noise σ without rheotaxis (F = 0). Each point is averaged over 20 runs. The total time of each simulation is tmax/(τ /I) = 10 4 with a time step dt/(τ /I) = 10 -2 .

FIG. 3 .

 3 FIG. 3. Time averaged orientations of each fish for different values of the noise without rheotaxis (F = 0). A short or large arrow indicates small and large polarization of a given fish respectively. a) For σ = 0.18, fish are well polarized. b) For σ = 0.27, a small zone of weak polarization appears. c) For a larger value σ = 0.54, only a few islands of polarized fish remain. d) Then, for σ = 1.0, no polarization remains (time averaged value ēi ≈ 0).

FIG. 4 .

 4 FIG. 4. Correlation function C as a function of the fish-to-fish distance d rescaled by R. The first part of the curves is associated with the correlation between close neighbor fish (d < 5a). Then, for d > 5a, an exponential decay is oberved. Orange circles: σ = 0.09; green left triangles: σ = 0.13 red squares: σ = 0.18; magenta down triangles: σ = 0.22; black diamonds: σ = 0.27; orange down triangles: σ = 0.36; green up triangles: σ = 0.45; purple crosses: σ = 0.54; brown plus: σ = 0.63; blue stars: σ = 0.72. Each solid line is a fit with an exponential decay such as exp(-d/λ), where λ is a correlation length (fitting parameter).

FIG. 6 .

 6 FIG. 6. Polarization P of the school against the flow as a function of F 1/2 /(στ 3/2 ) for different values of the parameter I/F . I/F = 20.0 (green crosses), I/F = 10.0 (magenta triangles), I/F = 2.0 (orange stars), I/F = 1.0 (blue diamond), I/F = 0.5 (red squares), I/F = 0 (black circles). A significant increase of P is observed when I/F increases. The arrow indicates the polarization difference between I/F = 20.0 and I/F = 0 at F 1/2 /(στ 3/2 ) = 1. Inset: Polarization difference of the school ∆P = P (I/F ) -P (I/F = 0). The colors and symbols represent the same I/F values as on the main figure. The presence of alignment interaction is very efficient for rheotaxis when I/F is large and F 1/2 /(στ 3/2 ) ≈ 1.

FIG. 8 .

 8 FIG. 8. Fluctuation of the school polarization as a function of polarization for different numbers of fish N . Green triangles N = 50; black diamonds N = 100; red squares N = 200; blue circles N = 400. Whatever N , the fluctuation maximum is around P = 0.6 ± 0.1.

2 and T 2 = 2 F 1 / 2 2

 22212 4 and < . > represents the averaging over the N fish. Let us call T 1 = στ 3/

ACKNOWLEDGEMENTS

We thank the Mission Interdisciplinaire du CNRS for a financial support. We thank our colleagues P. Moreau and P. Ballet for their technical help.

APPENDIX.

Different numbers of fish. We vary the number of fish, from N = 50 to N = 400 (see figure 7). Depending on N , the drop of polarization occurs in the range of noise 0.25 < σ < 0.5, however the amplitude of fluctuation always occurs around P = 0.6 ± 0.1 (Fig. 8). By changing N , we maintain the density ρ = N/(πR 2 ) constant and equal to ρ ≈ 0.5.

Dimensionless numbers

In the absence of rheotaxis (F = 0) the rescaling of eq.( 5) at stationnarity (see section C) leads to

where t = tI/τ and ω = ωτ /I. The dimensionless target is: