# Nonlinear Biomechanical Model of the Liver

Stéphanie Marchesseau, Simon Chatelin, Hervé Delingette

► To cite this version:

Stéphanie Marchesseau, Simon Chatelin, Hervé Delingette. Nonlinear Biomechanical Model of the Liver. Yohan Payan; Jacques Ohayon. Biomechanics of Living Organs. Hyperelastic Constitutive Laws for Finite Element Modeling, Elsevier, pp.243 - 265, 2017, 978-0-12-804009-6. 10.1016/b978-0-12-804009-6.00011-0. hal-03040194v2

HAL Id: hal-03040194

https://hal.science/hal-03040194v2

Submitted on 11 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non linear Biomechanical model of the Liver

Stéphanie Marchesseau<sup>a</sup>, Simon Chatelin<sup>b</sup>, Hervé Delingette<sup>c</sup>

<sup>a</sup>Clinical Imaging Research Centre, A*STAR-NUS, Singapore

<sup>b</sup>ICube, University of Strasbourg, CNRS, IHU Institute of Image Guided Surgery, Strasbourg, France

<sup>c</sup>Asclepios Research Project, INRIA Sophia Antipolis, France

## Abstract

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.

## 1. 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 [1]. 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 [2] 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 [3] 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 [4, 5], 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 [6], to perform augmented reality to visualize tumor locations during laparoscopic surgery [7, 8], to simulate the injection of gas inside the abdominal cavity (pneumoperitoneum) [9].
<figure><img src="image_2.png" /><figcaption>Figure 1: Liver Anatomy including the right and left lobes separated by the falciform ligament</figcaption></figure>

## 2. 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 venous 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.

## 3. 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.

### 3.1. 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 \cdot \frac{\partial W}{\partial 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, ($\lambda_i$, $I_i$ being the principal stretches and the invariants of the stretch tensor for $i = 1, 2, 3$ respectively).

<table><caption>Table 1: Main non-linear constitutive laws identified to model the liver tissue.</caption><tbody><tr><td rowspan="5">Polynomial forms</td><td>Neo-Hookean</td><td>$W = C_1(I_1 - 3)$</td></tr><tr><td>Mooney-Rivlin (generalized)</td><td>$W = \sum_{i+j>0}^{N} C_{ij}(I_1 - 3)^i(I_2 - 3)^j$</td></tr><tr><td>Ogden</td><td>$W = \sum_{k=1}^{N} \frac{\mu_k}{\alpha_k}(\lambda_1^{\alpha_k} + \lambda_2^{\alpha_k} + \lambda_3^{\alpha_k} - 3)$</td></tr><tr><td>Bogen</td><td>$W = \frac{\mu_1}{\alpha_1}(\lambda_1^{\alpha_1} + \lambda_2^{\alpha_1} + \lambda_3^{\alpha_1} - 1)$</td></tr><tr><td>Yeoh (reduced polynomial model)</td><td>$W = \sum_{k=1}^{N} C_k(I_1 - 3)^k$</td></tr><tr><td rowspan="3">Exponential and Logarithmic forms</td><td>Fung-Demiray</td><td>$W = \frac{C_1}{2C_2}(e^{C_2(I_1-3)} - 1)$</td></tr><tr><td>Logarithmic</td><td>$W = -C_1\ln(1 - C_2(\lambda_1^{\alpha_1} + \lambda_2^{\alpha_1} + \lambda_3^{\alpha_1} - 3))$</td></tr><tr><td>Exponential</td><td>$W = C_1(e^{C_2(\lambda_1^{\alpha_1} + \lambda_2^{\alpha_1} + \lambda_3^{\alpha_1})} - 1)$</td></tr><tr><td>Combined forms</td><td>Veronda-Westmann</td><td>$W = C_1(e^{C_3(I_1-3)} - 1) + C_2(I_2 - 3)$</td></tr></tbody></table>

#### 3.1.1. 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 [10, 11]. 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 ([12] (N = 1), [10] (N = 3), [13] (N = 2), [14] (N = 2), [15] (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 [10]. The use of high order Mooney-Rivlin constitutive laws is thought to result in the identification
of negative coefficients, which is not physically valid [16]. Reduced polynomial models [17], Bogen [10] and Ogden
strain functions ([16] (N = 3), [18] (N = 2 to 4), [17] (N = 4), [19] (N = 1), [20] (N = 1), [21](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.

3.1.2. Exponential and logarithmic hyperelastic models

Exponential and logarithmic constitutive forms have been presented as potential models representing the stress-
strain behavior of the hepatic tissue at least as well as the simple polynomial strain energy functions for both compres-
sion and elongation [10]. 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 re-
sponse of the hepatic tissue [22, 23, 24, 10, 25, 26]. Both exponential and logarithmic models show similar capacities
to simulate the hepatic tissue response [10].

3.1.3. Combined models

Combining the different forms of the strain energy functions is thought to result in the improvement of the model-
ing of the hepatic tissue, including high-order polynomial formulations [10]. The most popular approach consists in
combining polynomial with either logarithmic or exponential constitutive laws [14, 15], with the possibility of includ-
ing transverse isotropy using the fourth strain invariant [27]. 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 [28, 10].

3.1.4. 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 [29] and the first-order Mooney-Rivlin model [30, 31, 32]. This extension of polynomial
forms gives reasonable estimation of the hepatic tissue response at large strains over a large range of strain rates [30].
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 [33](N = 6). A non-linear strain-hardening fractional
derivative constitutive law has been proposed for liver tissue by Nicolle et al. [34]. 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⁻¹ to 0.7 s⁻¹). 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 [35]. While this Quasi-
Linear Viscoelastic (QLV) formulation has been successfully used to model the liver [23, 36, 25, 37, 38], 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 Kerdok in 2006 where she uses biphasic theory to include flow-independent hyperelasticity [39]. 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.

3.2. 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.

3.2.1. *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 experi-
ments. 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, vascu-
larization, 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) [20]. 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 [42, 39]. 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 [38, 33].

Up to now, the in vivo techniques used to identify the non-linear response of the liver tissue consisted either in
indentation [22, 42, 16, 31, 32, 38, 11, 17] or aspiration [36, 37] or Monitorized Endoscopic Grasper [26]. 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
[43]. The comparison between in vivo elastographic and ex vivo rheometric measurements has been investigated [44].
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.

3.2.2. 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⁻¹),
the effect of strain rate on porcine liver has been shown to be relatively insignificant by Chui et al. [10]. 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⁻¹), 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 [30]. 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⁻¹, respectively [20]. 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⁻¹ [19]. 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><img src="image_3.png" /><figcaption>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).</figcaption></figure>

### 3.2.3. 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 [27]. 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 [45]. Similarly, most of the experimental results and models in the literature are indicating isotropic behavior of the liver tissue.

### 3.2.4. 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 [29], swine [23, 24, 33] and bovine [45]). 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 Untaroiu et al. in 2013 and 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 [19, 21].

### 3.2.5. 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 [27], the capsule has specific characteristics, which are thought to result in modification of the global non-linear 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 [17]. 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 [46, 47]. 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 [48], confirming the importance of including the capsule to simulate the response of the entire liver.

### 3.3. 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 [49] discretized on linear tetrahedron finite elements which naturally leads to solving a linear system of equations whose inverse could eventually be precomputed [50].

Since linear elastic materials are not suitable for large displacements, several authors in the computer animation community have proposed corotational elastic models [51, 52] 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 [53] which exhibit a linear stress-strain relationship. Significant speed-up can be obtained by using reduced basis of deformation [54], Proper Generalized Decomposition [55] or by grouping expressions on edges, triangles and tetrahedra [56] 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 [57]. For instance, Miller et al. [58] 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) [59] 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.

## 4. 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 [41]. 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 hyperelasticity 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.

### 4.1. Optimised assembly of Finite Elements for Hyperelastic Materials

#### 4.1.1. 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
<figure><img src="image_4.png" /><figcaption>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$.</figcaption></figure>

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 [60], classical approaches involve the expensive computation of the second Piola-Kirchhoff (SPK) stress tensor S and the elasticity tensor $N_{ijkl} = 2\frac{\partial S_{ij}}{\partial 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 $\hat{N}$ can be fairly complex. Indeed, the first and second derivatives of the Jacobian $J = \det F$ 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, $\tilde{I}_1$ and $\tilde{I}_2$ particularly cumbersome and therefore computationally expensive to evaluate (see for instance [61] for details).

* The strain-displacement matrix $\hat{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 $\hat{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 \equiv 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.

### 4.1.2. 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ᵢ):

$$f_i = - \left( \frac{\partial W}{\partial Q_i} \right)^T \quad \text{and} \quad \mathbf{K}_{ij} = \left( \frac{\partial^2 W}{\partial Q_j \partial Q_i} \right)$$

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(\tilde{I})h^k(J)$$

or a function of it (for instance its exponential for Fung's law), where $\tilde{I} = (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) = \frac{dh^k}{dJ}$ requires a 1D derivation, and getting $S^k = 2\frac{\partial g^k(\tilde{I})}{\partial C}$ requires to combine well-known derivatives of the invariants (such as $\frac{\partial I_1}{\partial C} = \mathbf{I}d$ or $\frac{\partial I_2}{\partial C} = \mathbf{I}dI_1 - \mathbf{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 [41] for full derivation):

$$f_i = -V_0 \sum_k \left( h^{k'}(J) g^k(\tilde{I}) \left( \frac{\partial J}{\partial Q_i} \right)^T + h^k(J) \mathbf{F} S^k D_i \right)$$

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:

$$\frac{\partial^2 W(T_P)}{\partial Q_i \partial Q_j} = V_0 \sum_k (\mathbf{G}^k + \mathbf{H}^k + \mathbf{I}^k + \mathbf{\Lambda}^k + \mathbf{M}^k + \mathbf{R}^k)$$

where

$$\left\{ \begin{array}{lll} \mathbf{G}^k = \left( \frac{\partial h^{k'}(J)}{\partial Q_j} \right)^T \frac{\partial J}{\partial Q_i} g^k(\tilde{I}) & \mathbf{H}^k = h^{k'}(J) \frac{\partial}{\partial Q_j} \frac{\partial J}{\partial Q_i} g^k(\tilde{I}) & \mathbf{I}^k = h^{k'}(J) \left( \frac{\partial g^k(\tilde{I})}{\partial Q_j} \right)^T \frac{\partial J}{\partial Q_i} \\ \mathbf{\Lambda}^k = \left( \frac{\partial h^k(J)}{\partial Q_j} \right)^T D_i^T \mathbf{S}^k \mathbf{F}^T & \mathbf{M}^k = h^k(J) D_i^T \mathbf{S}^k \left( \frac{\partial \mathbf{F}^T}{\partial Q_j} \right)^T & \mathbf{R}^k = h^k(J) \left( \frac{\partial \mathbf{S}^k}{\partial Q_j} D_i \right)^T \mathbf{F}^T \end{array} \right.$$

We focus in more details on the term involving the third order tensor $\frac{\partial \mathbf{S}^k}{\partial Q_j}$:

$$\mathbf{R}^k = h^k(J) \left( \frac{\partial \mathbf{S}^k}{\partial Q_j} D_i \right)^T \mathbf{F}^T$$

The third order tensor $\frac{\partial \mathbf{S}^k}{\partial Q_j}$ is computed through the chain rule with the the elasticity tensor $\frac{\partial \mathbf{S}^k}{\partial C}$. Instead of representing the fourth order tensor $\frac{\partial \mathbf{S}^k}{\partial C}$ through 6×6 matrices, we define it through its application $\frac{\partial \mathbf{S}^k}{\partial C} : \mathbf{H}$ on any symmetric matrix $\mathbf{H}$. More precisely, we define the scalars $(\beta_l^k, \delta_l^k)$, and 3 × 3 symmetric matrices $(\mathbf{A}_l^k, \mathbf{B}_l^k)$ such that :

$$\frac{\partial \mathbf{S}^k}{\partial \mathbf{C}} : \mathbf{H} = \sum_l \beta_l^k \mathbf{A}_l^k \mathbf{H} \mathbf{A}_l^k + \delta_l^k (\mathbf{H} : \mathbf{B}_l^k) \mathbf{B}_l^k$$
where $\mathbf{A}: \mathbf{B} = \text{tr}(\mathbf{B}^T \mathbf{A})$ for any two matrices $\mathbf{A}, \mathbf{B}$. Finally, the term $\mathbf{R}^k$ is found to be a combination of two terms:

$$h^k(J)\mathbf{F}\mathbf{L}_i^k(i, j)\mathbf{F}^T \text{ and } h^k(J)\mathbf{F}\mathbf{U}_i^k(i, j)\mathbf{F}^T$$

where $\mathbf{L}_i^k(i, j)$ and $\mathbf{U}_i^k(i, j)$ are linear matrices:

$$\left\{ \begin{aligned} \mathbf{L}_i^k(i, j) &= \beta_i^k (\mathbf{A}_i^k D_i \otimes D_j \mathbf{A}_i^k + \mathbf{A}_i^k (D_j \cdot \mathbf{A}_i^k D_i)) \\ \mathbf{U}_i^k(i, j) &= \delta_i^k (\mathbf{B}_i^k D_j \otimes D_i \mathbf{B}_i^k) \end{aligned} \right.$$

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 $\mathbf{A}_i^k$ and $\mathbf{B}_i^k$ are constant and therefore matrices $\mathbf{L}_i^k(i, j)$ and $\mathbf{U}_i^k(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, [62] perform an SVD decomposition of the deformation gradient matrix. To avoid this computationally expensive decomposition, we propose instead to regularize the term $\mathbf{G}^k = h^{k''}(J) g^k(\tilde{I}) \frac{\partial J}{\partial Q_j} \otimes \frac{\partial J}{\partial Q_i}$ by replacing it with the following expression:

$$\mathbf{G}^k = h^{k''}(J)g^k(\tilde{I}) \left( (1 - \epsilon) \frac{\partial J}{\partial Q_j} \otimes \frac{\partial J}{\partial Q_i} + \epsilon \frac{\partial J}{\partial Q_j} \cdot \frac{\partial J}{\partial Q_i} \mathbf{Id} \right)$$

The closer $\epsilon$ is to 1 the closer the $\mathbf{G}^k$ matrix is to a diagonal matrix. In practice, we set $\epsilon = (1 - J)$ if $0 \le J \le 1$, $\epsilon = 0$ if $J \ge 1$ and $\epsilon = 1$ if $J \le 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 it does not directly compute the (first or) second Piola Kirchhoff stress tensor $\mathbf{S} = 2\frac{W}{\partial C}$. If such stress tensors are required, then they need to be computed at extra cost (see below).

### 4.2. Visco-hyperelasticity based on Prony series

To model accurately the viscoelasticity of the liver, we propose to rely on Prony series [59]. This consists in adding to the hyperelastic SPK stress tensor now denoted $\mathbf{h}\mathbf{S}$ some time dependent stresses. This time dependence is given by $\alpha(t) = \alpha_\infty + \sum_i \alpha_i e^{-t/\tau_i}$ with the condition $(\alpha_\infty + \sum_i \alpha_i) = 1$. The visco-hyperelastic SPK tensor noted $\mathbf{v}\mathbf{S}$ can then be written as:

$$\mathbf{v}\mathbf{S} = \int_0^t \alpha(t-t') \frac{\partial \mathbf{h}\mathbf{S}}{\partial t'} dt' = \mathbf{h}\mathbf{S} - \sum_i \gamma_i \quad \text{where} \quad \gamma_i = \int_0^t \alpha_i (1 - e^{(t'-t)/\tau_i}) \frac{\partial \mathbf{h}\mathbf{S}}{\partial t'} dt'$$

After a discretization over time this results in the recursive formula between time $n-1$ and time $n$:

$$\gamma_i^n = a_i \mathbf{h}\mathbf{S}^n + b_i \gamma_i^{n-1} \quad \text{where} \quad a_i = \frac{\Delta t \alpha_i}{\Delta t + \tau_i} \quad \text{and} \quad b_i = \frac{\tau_i}{\Delta t + \tau_i}$$

$\Delta 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 $\mathbf{h}\mathbf{S}^n$. This is done for each time step $n$ (dropping the index for clarity) by computing the inverse deformation gradient:

$$\mathbf{h}\mathbf{S} = \mathbf{F}^{-1} \left( \sum_k (f^{k'}(J)g^k(\tilde{I})\mathbf{J}\mathbf{F}^{-T} + f^k(J)\mathbf{F}\mathbf{h}\mathbf{S}^k) \right) \quad \text{where} \quad \mathbf{F}^{-1} = \left( \sum_{l=1}^4 P_l \otimes \frac{\partial J}{\partial Q_l} \right) / J$$

The visco-hyperelastic nodal forces $\mathbf{v}f$ are therefore related to the hyperelastic ones $\mathbf{h}f$ by

$$\mathbf{v}f_i = \mathbf{h}f_i + V_0 \mathbf{F} \sum_i \gamma_i D_i$$
Moreover, once we have $\gamma_i^{n-1}$ the stiffness matrix is also slightly updated at each time $n$ from its hyperelastic formulation $^hK$:

$$ ^vK_{ij}^n = ^hK_{ij}^n \left( 1 - \sum_k a_k \right) - V_0 D_j^T \left( \sum_k b_k \gamma_k^{n-1} \right) D_i \mathbf{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 $\gamma_i^n$ and $F^{-1}$.

### 4.3. Poro-elasticity

We follow Kerdok's porosity model [39] and consider the liver as a fluid-filled sponge. This model is made of a volumetric component (represented by $\sigma_{Heq}$ Figure 3) and a fluid phase with variable $P_{fluid}$. The volumetric component is governed by Hencky's elasticity [63] 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:

$$ \sigma_{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 $1 - r_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 $\sigma_{Heq}$ when $J \le J_0$ by its tangent curve at $J_0$ (see Figure 3(Right)). We set $J_0 = 1 - r_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_{fluid}$ is governed by the variation of volume change and a diffusive process:

$$ \frac{1}{K_{lim}} \dot{P}_{fluid} = \kappa \nabla^2 P_{fluid} - \frac{\dot{J}}{J} $$

where $\kappa$ is the permeability parameter. In Kerdok's model, the permeability $\kappa$ 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:

$$ \sigma_p = \sigma_{heq} \mathbf{Id} - P_{fluid} \mathbf{Id} $$

The Cauchy stress is translated as a poro-elastic force:

$$ {}^p f_i = -\sigma_p \left( \frac{\partial J}{\partial Q_i} \right)^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_{ij} = V_0 \left( \frac{\partial \sigma_p}{\partial J} \left( \frac{\partial J}{\partial Q_j} \right)^T \left( \frac{\partial J}{\partial Q_i} \right) + \sigma_p \frac{\partial^2 J}{\partial Q_j \partial Q_i} \right) $$

## 5. Results

### 5.1. MJED validation

First, we compared the computation time between the MJED and the classical FEM method, referred to as "Standard FEM", implemented in SOFA¹. 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 [41]. We simulated the deformation of
<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of the computation times of nodal forces and stiffness matrices between two different discretization methods averaged over 100 iterations.</figcaption></figure>

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)<fnref n="2" /> where several elastic and hyperelastic materials are implemented and against analytical solutions for simple cube model. The mean relative difference is around 10⁻⁶ 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.

## 5.2. 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^*(\omega) = G'(\omega) + iG''(\omega)$, 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_{\infty} + g_1e^{-t/\tau_1} + g_2e^{-t/\tau_2})$$

<footnote n="1">¹SOFA is an Open Source medical simulation software available at www.sofa-framework.org</footnote>

<footnote n="2">²FEBio is an opensource software package for FE analysis available at mrl.sci.utah.edu</footnote>
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. [17, 23, 31, 11, 15]

5.3. *Model Parameter fitting from experiments*

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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⁻¹) [17, 23, 31, 11, 15]. 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.
<table><caption>Table 2: Values of the parameters identified from the DMA Tests for the viscosity and from the literature for the porosity.</caption><thead><tr><th colspan="5">Viscosity</th><th colspan="4">Porosity</th></tr><tr><th>G<sub>0</sub>(Pa)</th><th>g<sub>1</sub>(Pa)</th><th>τ<sub>1</sub>(s)</th><th>g<sub>2</sub>(Pa)</th><th>τ<sub>2</sub>(s)</th><th>r<sub>w</sub></th><th>K<sub>0</sub>(Pa)</th><th>K<sub>lim</sub>(kPa)</th><th>κ(m<sup>4</sup>/Ns)</th></tr></thead><tbody><tr><td>770</td><td>0.235</td><td>0.27</td><td>0.333</td><td>0.03</td><td>0.5</td><td>400</td><td>2.2</td><td>20</td></tr></tbody></table>

## 5.4. 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 $1 - \sum a_k$ at infinite 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 [39] 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 $\dot{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.

# 6. 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 visco-poro-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 [64] due to the excessive accumulation of collagen fibers. Fatty liver disease (hepatic steatosis) also translates into a significant increase in liver stiffness [67]. Elastography imaging techniques based on Ultrasound [65, 66] or Magnetic Resonance Imaging [68] 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
<figure><img src="image_7.png" /><figcaption>Figure 6: (Top) Pressure field of the porous component on a liver under gravity, during one oscillation. (Bottom) Addition of viscosity to hyper-elasticity: Comparison of the maximum amplitudes and final states. (Black) Initial position, (Blue) Hyperelastic liver, (Pink) Visco-hyperelastic liver.</figcaption></figure>

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.

References

[1] C. D. Untaroiu, Y.-C. Lu, S. K. Siripurapu, A. R. Kemper, Modeling the biomechanical and injury response of human liver parenchyma under tensile loading, Journal of the Mechanical Behavior of Biomedical Materials 41 (2015) 280–291.

[2] F. Nickel, J. A. Brzoska, M. Gondan, H. M. Rangnick, J. Chu, H. G. Kenngott, G. R. Linke, M. Kadmon, L. Fischer, B. P. Muller-Stich, Virtual reality training versus blended learning of laparoscopic cholecystectomy: a randomized controlled trial with laparoscopic novices, Medicine (Baltimore) 94 (20) (2015) e764.

[3] H. Delingette, N. Ayache, Soft tissue modeling for surgery simulation, in: N. Ayache (Ed.), Computational Models for the Human Body, Handbook of Numerical Analysis (Ed : Ph. Ciarlet), Elsevier, 2004, pp. 453–550.

[4] O. Oktay, L. Zhang, T. Mansi, P. Mountney, P. W. Mewes, S. Nicolau, L. Soler, C. Chefd'Hotel, Biomechanically driven registration of pre-to intra-operative 3d images for laparoscopic surgery, in: Medical Image Computing and Computer-Assisted Intervention - MICCAI 2013 - 16th International Conference, Nagoya, Japan, September 22-26, 2013, Proceedings, Part II, 2013, pp. 1-9.

[5] A. L. Simpson, P. Dumpuri, W. R. Jarnagin, M. I. Miga, Soft Tissue Biomechanical Modeling for Computer Assisted Surgery, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, Ch. Model-Assisted Image-Guided Liver Surgery Using Sparse Intraoperative Data, pp. 7–40.

[6] Y. Kobayashi, H. Watanabe, T. Hoshi, K. Kawamura, M. G. Fujie, Soft Tissue Biomechanical Modeling for Computer Assisted Surgery, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, Ch. Viscoelastic and Nonlinear Liver Modeling for Needle Insertion Simulation, pp. 41–67.

[7] R. Plantefève, I. Peterlik, N. Haouchine, S. Cotin, Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery, Annals of Biomedical Engineering 44 (1) (2015) 139-153.
[8] N. Haouchine, S. Cotin, I. Peterlík, J. Dequidt, M. Sanz-Lopez, E. Kerrien, M. Berger, Impact of soft tissue heterogeneity on augmented reality for liver surgery, IEEE Trans. Vis. Comput. Graph. 21 (5) (2015) 584-597.

[9] J. Bano, A. Hostettler, S. Nicolau, S. Cotin, C. Doignon, H. S. Wu, M. H. Huang, L. Soler, J. Marescaux, Simulation of pneumoperitoneum for laparoscopic surgery planning, in: Medical Image Computing and Computer-Assisted Intervention - MICCAI 2012 - 15th International Conference, Nice, France, October 1-5, 2012, Proceedings, Part I, 2012, pp. 91-98.

[10] C. Chui, E. Kobayashi, X. Chen, T. Hisada, I. Sakuma, Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation, Medical and Biological Engineering and Computing 42 (6) (2004) 787-798.

[11] B. Ahn, J. Kim, Measurement and characterization of soft tissue behavior with surface deformation and force response under large deformations, Medical Image Analysis 14 (2) (2010) 138-148.

[12] T. Hu, J. P. Desai, A biomechanical model of the liver for reality-based haptic feedback, in: Medical Image Computing and Computer-Assisted Intervention-MICCAI 2003, Springer, 2003, pp. 75-82.

[13] S. Umale, C. Deck, N. Bourdet, P. Dhumane, L. Soler, J. Marescaux, R. Willinger, Experimental mechanical characterization of abdominal organs: liver, kidney & spleen, Journal of the mechanical behavior of biomedical materials 17 (2013) 22-33.

[14] Y. Fu, C. Chui, C. Teo, Liver tissue characterization from uniaxial stress-strain data using probabilistic and inverse finite element methods, Journal of the mechanical behavior of biomedical materials 20 (2013) 105-112.

[15] Y. Fu, C. Chui, Modelling and simulation of porcine liver tissue indentation using finite element method and uniaxial stress-strain data, Journal of biomechanics 47 (10) (2014) 2430-2435.

[16] T. Hu, J. P. Desai, Characterization of soft-tissue material properties: large deformation analysis, in: Medical Simulation, Springer, 2004, pp. 28-37.

[17] K. Lister, Z. Gao, J. P. Desai, Development of in vivo constitutive models for liver: application to surgical simulation, Annals of biomedical engineering 39 (3) (2011) 1060-1073.

[18] Z. Gao, J. P. Desai, Estimating zero-strain states of very soft tissue under gravity loading using digital image correlation, Medical image analysis 14 (2) (2010) 126-137.

[19] C. D. Untaroiu, Y.-C. Lu, Material characterization of liver parenchyma using specimen-specific finite element models, Journal of the mechanical behavior of biomedical materials 26 (2013) 11-22.

[20] Y.-C. Lu, A. R. Kemper, C. D. Untaroiu, Effect of storage on tensile material properties of bovine liver, Journal of the mechanical behavior of biomedical materials 29 (2014) 339-349.

[21] C. D. Untaroiu, Y.-C. Lu, S. K. Siripurapu, A. R. Kemper, Modeling the biomechanical and injury response of human liver parenchyma under tensile loading, Journal of the mechanical behavior of biomedical materials 41 (2015) 280-291.

[22] F. J. Carter, T. G. Frank, P. J. Davies, D. McLean, A. Cuschieri, Measurements and modelling of the compliance of human and porcine organs, Medical Image Analysis 5 (4) (2001) 231-236.

[23] A. Tamura, K. Omori, K. Miki, J. B. Lee, K. H. Yang, A. I. King, Mechanical characterization of porcine abdominal organs., Stapp car crash journal 46 (2002) 55-69.

[24] J. D. Brown, J. Rosen, Y. S. Kim, L. Chang, M. N. Sinanan, B. Hannaford, In-vivo and in-situ compressive properties of porcine abdominal soft tissues, Studies in Health Technology and Informatics (2003) 26-32.

[25] E. Roan, K. Vemaganti, The nonlinear material properties of liver tissue determined from no-slip uniaxial compression experiments, Journal of biomechanical engineering 129 (3) (2007) 450-456.

[26] J. Rosen, J. D. Brown, S. De, M. Sinanan, B. Hannaford, Biomechanical properties of abdominal organs in vivo and postmortem under compression loads, Journal of Biomechanical Engineering 130 (2) (2008) 021020.

[27] C. Chui, E. Kobayashi, X. Chen, T. Hisada, I. Sakuma, Transversely isotropic properties of porcine liver tissue: experiments and constitutive modelling, Medical & biological engineering & computing 45 (1) (2007) 99-106.

[28] H. Yin, L. Sun, G. Wang, M. W. Vannier, Modeling of elastic modulus evolution of cirrhotic human liver, Biomedical Engineering, IEEE Transactions on 51 (10) (2004) 1854-1857.

[29] K. Miller, Constitutive modelling of abdominal organs, Journal of biomechanics 33 (3) (2000) 367-373.

[30] Z. Liu, L. Bilston, Large deformations shear properties of liver tissue, Biorheology 39 (6) (2002) 735-742.

[31] J. Kim, M. A. Srinivasan, Characterization of viscoelastic soft tissue properties from in vivo animal experiments and inverse fe parameter estimation, in: Medical Image Computing and Computer-Assisted Intervention-MICCAI 2005, Springer, 2005, pp. 599-606.

[32] E. Samur, M. Sedef, C. Basdogan, L. Avtan, O. Duzgun, A robotic indenter for minimally invasive measurement and characterization of soft tissue response, Medical Image Analysis 11 (4) (2007) 361-373.

[33] F. Sato, Y. Yamamoto, D. Ito, J. Antona-Makoshi, S. Ejima, K. Kamiji, T. Yasuki, Hyper-viscoelastic response of perfused liver under dynamic compression and estimation of tissue strain thresholds with a liver finite element model, in: IRCOBI conference, 2013.

[34] S. Nicolle, P. Vezin, J.-F. Palierne, A strain-hardening bi-power law for the nonlinear behaviour of biological soft tissues, Journal of biomechanics 43 (5) (2010) 927-932.

[35] Y. Fung, Elasticity of soft tissues in simple elongation, American Journal of Physiology-Legacy Content 213 (6) (1967) 1532-1544.

[36] A. Nava, E. Mazza, O. Haefner, M. Bajka, Experimental observation and modelling of preconditioning in soft biological tissues, in: Medical Simulation, Springer, 2004, pp. 1-8.

[37] A. Nava, E. Mazza, M. Furrer, P. Villiger, W. Reinhart, In vivo mechanical characterization of human liver, Medical image analysis 12 (2) (2008) 203-216.

[38] P. Jordan, S. Socrate, T. Zickler, R. Howe, Constitutive modeling of porcine liver in indentation using 3d ultrasound imaging, Journal of the mechanical behavior of biomedical materials 2 (2) (2009) 192-201.

[39] A. E. Kerdok, Characterizing the nonlinear mechanical response of liver to surgical manipulation, Ph.D. thesis, Harvard University (2006).

[40] S. Raghunathan, D. Evans, J. L. Sparks, Poroviscoelastic modeling of liver biomechanical response in unconfined compression, Annals of biomedical engineering 38 (5) (2010) 1789-1800.

[41] S. Marchesseau, T. Heimann, S. Chatelin, R. Willinger, H. Delingette, Fast porous visco-hyperelastic soft tissue model for surgery simulation: Application to liver surgery, Progress in biophysics and molecular biology 103 (2) (2010) 185-196.
[42] M. Ottensmeyer, Tempest 1-d: an instrument for measuring solid organ soft tissue properties, Experimental Techniques 26 (3) (2002) 48–50.

[43] J. Vappou, Magnetic resonance- and ultrasound imaging- based elasticity imaging methods: A review, Critical Reviews in Biomedical Engineering 40 (2).

[44] S. Chatelin, J. Oudry, N. Périchon, L. Sandrin, P. Allemann, L. Soler, R. Willinger, In vivo liver tissue mechanical properties by transient elastography: Comparison with dynamic mechanical analysis, Biorheology 48 (2) (2011) 75–88.

[45] F. Pervin, W. W. Chen, T. Weerasooriya, Dynamic compressive response of bovine liver tissues, Journal of the mechanical behavior of biomedical materials 4 (1) (2011) 76–84.

[46] A. Brunon, K. Bruyere-Garnier, M. Coret, Mechanical characterization of liver capsule through uniaxial quasi-static tensile tests until failure, Journal of biomechanics 43 (11) (2010) 2221–2227.

[47] A. Brunon, K. Bruyere-Garnier, M. Coret, Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests, Journal of the mechanical behavior of biomedical materials 4 (8) (2011) 1572–1581.

[48] S. Umale, S. Chatelin, N. Bourdet, C. Deck, M. Diana, P. Dhumane, L. Soler, J. Marescaux, R. Willinger, Experimental in vitro mechanical characterization of porcine glisson's capsule and hepatic veins, Journal of biomechanics 44 (9) (2011) 1678–1683.

[49] S. Cotin, H. Delingette, N. Ayache, A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation, The Visual Computer 16 (8) (2000) 437–452.

[50] D. L. James, D. K. Pai, ArtDefo accurate real time deformable objects, in: Computer Graphics (SIGGRAPH), 1999, pp. 65–72.

[51] S. Suwelack, S. Röhl, R. Dillmann, A.-L. Wekerle, H. Kenngott, B. Müller-Stich, C. Alt, S. Speidel, Computational Biomechanics for Medicine: Deformation and Flow, Springer New York, New York, NY, 2012, Ch. Quadratic Corotated Finite Elements for Real-Time Soft Tissue Registration, pp. 39–50.

[52] H. Courtecuisse, H. Jung, J. Allard, C. Duriez, D. Y. Lee, S. Cotin, Gpu-based real-time soft tissue deformation with cutting and haptic feedback, Progress in Biophysics and Molecular Biology 103 (23) (2010) 159 – 168, special Issue on Biomechanical Modelling of Soft Tissue Motion.

[53] H. Delingette, N. Ayache, Soft tissue modeling for surgery simulation, in: Computational Models for the Human Body, Elsevier, 2004, pp. 453–550.

[54] J. Barbič, D. L. James, Real-time subspace integration for St. Venant-Kirchhoff deformable models, ACM TOG (SIGGRAPH 2005) 24 (3) (2005) 982–990.

[55] S. Niroomandi, D. Gonzlez, I. Alfaro, F. Bordeu, A. Leygue, E. Cueto, F. Chinesta, Real-time simulation of biological soft tissues: a pgd approach, International Journal for Numerical Methods in Biomedical Engineering 29 (5) (2013) 586–600.

[56] G. Picinbono, H. Delingette, N. Ayache, Non-Linear Anisotropic Elasticity for Real-Time Surgery Simulation, Graphical Models 65 (5) (2003) 305–321.

[57] S. Speidel, S. Roehl, S. Suwelack, R. Dillmann, H. Kenngott, B. Mueller-Stich, Intraoperative surface reconstruction and biomechanical modeling for soft tissue registration, in: Proc. Joint Workshop on New Technologies for Computer/Robot Assisted Surgery, 2011, pp. 0–0.

[58] K. Miller, G. Joldes, D. Lance, A. Wittek, Total lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation, Communications in Numerical Methods in Engineering 23 (2006) 121–134.

[59] Z. Taylor, O. Comas, M. Cheng, J. passenger, D. hawkes, D. . Atkinson, S. Ourselin, On modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and gpu execution, medical image Analysis 13 (2009) 234–244.

[60] C. Zienkiewicz, R. Taylor, The Finite Element Method, Volume 2 : Solid Mechanics, Butterworth-Heinemann, 2000.

[61] J. Weiss, B. Maker, S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Computer methods in applied mechanics and engineering 135 (1996) 107–128.

[62] J. Teran, E. Sifakis, S. S. Blemker, V. Ng-Thow-Hing, C. Lau, R. Fedkiw, Creating and simulating skeletal muscle from the visible human data set, IEEE Transactions on Visualization and Computer Graphics 11 (3) (2005) 317–328.

[63] H. Xiao, L. Chen, Hencky's elasticity model and linear stress-strain relations in isotropic finite hyperelasticity, Acta Mechanica 157 (2002) 51–60.

[64] L. Castra, J. Vergniol, J. Foucher, B. L. Bail, E. Chanteloup, M. Haaser, M. Darriet, P. Couzigou, V. de Ldinghen, Prospective comparison of transient elastography, fibrotest, apri, and liver biopsy for the assessment of fibrosis in chronic hepatitis c, Gastroenterology 128 (2) (2005) 343 – 350.

[65] J.-L. Gennisson, T. Deffieux, M. Fink, M. Tanter, Ultrasound elastography: principles and techniques, Diagnostic and interventional imaging 94 (5) (2013) 487–495.

[66] M. Fatemi, J. F. Greenleaf, Ultrasound-stimulated vibro-acoustic spectrography, Science 280 (5360) (1998) 82–85.

[67] L. Huwart, C. Sempoux, E. Vicaut, N. Salameh, L. Annet, E. Danse, F. Peeters, L. C. ter Beek, J. Rahier, R. Sinkus, et al., Magnetic resonance elastography for the noninvasive staging of liver fibrosis, Gastroenterology 135 (1) (2008) 32–40.

[68] R. Muthupillai, D. Lomas, P. Rossman, J. Greenleaf, A. Manduca, R. Ehman, Magnetic resonance elastography by direct visualization of propagating acoustic strain waves, Science 269 (5232) (1995) 1854–1857.