Parameters identification method for breast biomechanical numerical model

by   Diogo Lopes, et al.

Bio-mechanical breast simulations are based on a gravity free geometry as a reference domain and a nonlinear mechanical model parameterised by physical coefficients. As opposed to complex models proposed in the literature based on medical imagery, we propose a simple but yet realistic model that uses a basic set of measurements easy to realise in the context of routinely operations. Both the mechanical system and the geometry are controlled with parameters we shall identify in an optimisation procedure. We give a detailed presentation of the model together with the optimisation method and the associated discretisation. Sensitivity analysis is then carried out to evaluate the robustness of the method.


Discussion on Mechanical Learning and Learning Machine

Mechanical learning is a computing system that is based on a set of simp...

Inverse analysis of material parameters in coupled multi-physics biofilm models

In this article we propose an inverse analysis algorithm to find the bes...

An Hybrid Method for the Estimation of the Breast Mechanical Parameters

There are several numerical models that describe real phenomena being us...

Spontaneous Fruit Fly Optimisation for truss weight minimisation: Performance evaluation based on the no free lunch theorem

Over the past decade, several researchers have presented various optimis...

Identifying Mechanical Models through Differentiable Simulations

This paper proposes a new method for manipulating unknown objects throug...

A relativistic extension of Hopfield neural networks via the mechanical analogy

We propose a modification of the cost function of the Hopfield model who...

An electro-chemo-mechanical framework for predicting hydrogen uptake in metals due to aqueous electrolytes

We present a theoretical and numerical scheme that enables quantifying h...

1 Introduction

Bio-mechanical breast modelling and simulation with the finite element method (del Palomar et al., 2008, Eder et al., 2014) provides powerful analysis tools to predict 3-D breast deformations for various bio-medical applications. Breast tissue modelling is usually based on non-linear mechanical models such as the Neo-hookean model expressed in the stress-free unloaded reference state of the breast.

Both the parameters of the mechanical model and the initial gravity-free configuration are unknown (Rajagopal et al., 2007) and have to be identified for each patient. To solve the difficulty in determining the material parameters of soft tissues for simulations and providing the reference configuration, optimisation procedures are usually carried out on the basis of MR image or radiological acquisitions such as mammographic plate compression (Lorenzen et al., 2003, Eder et al., 2014). Breast models have also been assessed using predicted location of anatomical landmarks and selected in breast images acquired before and after in vivo compression by visual comparison. Then a complex procedure of reverse problem is achieved (Affagard et al., 2015, Isvilanonda et al., 2016, Eder et al., 2014) for determining the reference state and the associated coefficients.

The methods proposed in the literature (del Palomar et al., 2008, Eder et al., 2014) require images acquisition, complex identification procedures with high computational cost (Han et al., 2012, Eiben et al., 2016) that would be unnecessary or too expensive for routinely or low-cost practices. From a practical point of view, quality does not necessarily mean high accuracy or high amounts of detailed information like the models proposed in (Azar et al., 2002), but to manage relevant and satisfactory results for surgeons. In this context, we propose an alternative and simpler method that provides an effective numerical technology enable to accurately predict breast deformations and running on a mid range laptop within few minutes following the model proposed by Cardoso et al. (2003). The key ideas are, on the one hand, the introduction of a simple set of in vivo breast measures easy to achieve, and on the other hand, a four-parameters non-linear mechanical model defining on a simplified two-parameters reference configuration (the one without gravity). Measurements provide the data we use to fit the model parameters (both mechanical and geometrical). Given an a priori

geometry with only two degrees of freedom reduces the inverse problem difficulty,

i.e. the recovery of the stress-free unloaded reference state and provides a faster identification procedure.

Concerning related works, in Cardoso et al. (2003), the authors treat a similar model but only consider a two-parameter problem, namely the elasticity coefficients of the core domain, while the present contribution treats both the geometrical and mechanical aspects including the skin effect leading to a more complete six parameters model. The use of nonlinear models to simulate the skin, the muscles, and the tissues, is well-developed in the bio-mechanics context (Gascón et al., 2014, Fung, 2012, Rajagopal et al., 2004) while the finite element method is a popular technique for the discretisation (Payan, 2012, Zhong et al., 2005, Lee and Zhuang, 2012) using alternatively the weak formulation or the minimisation framework (Ball, 1983, Ciarlet, 1998, Picinbono et al., 2003). The Neo-hookean is considered as a well-adapted mechanical representation of the breast (del Palomar et al., 2008) and more generally for human tissues (Affagard et al., 2015, Isvilanonda et al., 2016). For the particular case of breast, we differentiate the skin from the core assuming a couple of coefficients for each material.

The paper is organised as follows. We present the breast model with its mathematical formulation in Section 2, followed by the parameter identification method in Section 3. In Section 4 we show the discretisation of this model and parameter evaluation method while results obtained using synthetic measurements are presented in Section 5. Finally, in Section 6 conclusions are presented.

2 Breast Modelling

The breast is a complex structure constituted of a mass of glandular tissue encased in variable quantities of fat, that account for its characteristic round shape, connected to the skin trough a series of ligaments. All these types of tissues possess different mechanical properties while the proportion of these materials varies with factors like genetics and age (Poplack et al., 2004). Moreover, the breast skin plays a major role as a wrapper tissue enjoying specific mechanical properties.

Deformation evaluation of the breast over external actions is achieved by considering a simplified stress-free geometrical domain of the breast equipped with the Neo-hookean mechanical model (del Palomar et al., 2008) where we differentiate the breast core and skin. Additionally, we also introduce the Chassaignac space (Cardoso et al., 2010) to reproduce the breast mobility. The gravity-free reference breast domain is a piece of spherical cap where the plane section is attached to the torso. The domain is parameterised with the radius of the sphere while represents the non-truncated length as displayed in Figure 1. The cap which is placed in the torso plane corresponds to a circle of radius which is smaller than .

Figure 1: Breast geometry configurations. Radius characterised the main sphere while identifies the truncated part following the negative direction.

The final and more realistic shape of the breast is obtained by determining the effect of gravity on the breast tissues as observed in figure 2.

To prescribe the boundary condition and compute the energy associated to the skin, we introduce the following notations, reproduced on Figure 2: is the back side plane of the breast which attaches to the torso, represents the surface associated to the skin of the breast and is the arc of the infra-mammary fold.

Figure 2: Notation and geometry of the gravity-free geometry (left) and configuration in the gravity field (right)

The hyperelastic neo-Hookean compressible relations (Cardoso et al., 2010) equipped with the so-called Chassaignac space represent a well-accepted model for biomechanical soft tissues. The Chassaignac space corresponds to a mobile zone located between the breast and the trunk which acts as a spring to maintain the breast close to the trunk. Under the gravity, the breast corresponds to a minimisation of the energy functional (Ball, 1983, Ciarlet, 1998, Picinbono et al., 2003) associated with the neo-Hookean system. Such a functional aggregates the internal energy due to the bulk, the energy deriving from the skin displacement, the gravitational energy characterised by the gravitational field and the spring energy associated to Chassaignac space. Anisotropic model for skin Groves et al. (2013) is also investigated to provide a more sophisticated description taking into account the mechanical behaviour of skin under large deformations.

For any generic point

, we seek the new position vector field

which minimises the energy functional given by


subject to the constraints


where is the Chassaignac coefficient, and stands for the bulk density. Constraint (2) represents the fixation of the infra-mammary fold on the torso, i.e., any point in maintain their position while constraint (3) states that any point on the breast plane attached to the torso (Chassaignac space) only moves inside that trunk plane.

The expressions for the volume strain-energy and the skin strain-energy densities, respectively represented by and , are given by

where is the Jacobian matrix of and are the Lamé parameters for the breast, and

where are the Lamé parameters for the skin while is the Jacobi matrix of the superficial (skin) displacement ( is the restriction of on using local two-parameters representation since is a surface).

3 Parameters identification

Let denote the six parameters which represents the geometrical and the physical degrees of freedom respectively. First, for the two geometrical parameters , we create the stress-free configuration , as shown in Figure 1, which defines the operator


Secondly, for a set of physical parameters , minimisation of the energy over the set of regular functions defined on provides the solution . At last, applying the displacement over the reference domain provides the deformed domain which, at the end of the day, defines the operator


Third, to perform the parameters identification, we consider the following measurements: the volume, the surface area of the breast’s skin and three sizes, namely , and , as displayed in Figure 3.

Figure 3: Size characteristic (, and ) of a breast (stand-up position)

Measures are performed for three patient positions, stand up, inclined (45 degrees) and horizontal, providing 15 real values denoted . On the other hand, given a set of parameters and solving three times the direct problem corresponding to the three positions provides the vectors . The problem consists in seeking the set of parameters such that is equal to in the least square sense. Introducing the cost function


we seek for which minimises the errors between the experimental data and the theoretical solution values.

To provide the optimal set of parameters, an iterative process is considered by evaluating the sequences such that . We use the Gauss-Newton method to take advantage that function casts in the specific form


The first-order Taylor series expansion reads


where is the Jacobian matrix. From relations 7 and 8 we deduce the first term of error for a perturbation of the parameters

Thus, by identification, we deduce


We seek a perturbation that minimises , i.e., . Therefore Equation 9 provides the expression


Since a direct computation of the jacobian matrix is not possible, we numerically evaluate each partial derivative setting

where is fixed by the user in function of the minimisation problem and .

4 Discretisation

To carry out numerical simulations, we introduce the discretisation of the breast model and the optimisation problem. We denote by a mesh of the gravity free domain constituted of non-overlapping tetrahedron cells , , and vertices , while

stands for the discrete domain. Moreover, , , represents the faces of the tetrahedrons of the mesh that belong to . Quantities and represent the volume and the area of the cell and the triangle respectively. We also use a local indexation and denote by , , the vertices of and by , , the vertices of . To discretise function , we associate to each node an approximation and denote by the continuous, linear piecewise function while vector collects the components of the new positions.

Vector corresponds to the new configuration but not all the entries are necessarily unknowns of the problem since some of them are characterised by the boundary conditions. The sub-vector of only contains the unknown values we shall use in the minimisation process while the boundary conditions define an operator

which provides the other entries to complete vector . In the present contribution, condition (2) yields that for any while relation (3) implies that for any , we set to maintain interface on the trunk plane . This two conditions completely define the vector of unknowns and operator .

4.1 The energy functional

Discrete version of the energy functional is given by

where , , and represent the volume energy, the surface energy, the energy from the gravitational displacement, and the energy associated to the Chassaignac space respectively.

For a new configuration characterised by the approximation and stored in vector , the internal energy on tetrahedron is given by

where is the matrix solution of the linear system

finally having .

For the surface energy, the discrete piecewise linear function transforms a triangle with vertices into a triangle with vertices . The method is similar to the surface variation using in Lee et al. (2018) to evaluate the energy deriving from the displacement variations. Since the translation and the rotation do not change the stress due to the deformation, we assume that is collinear to and belongs to the same plane as triangle . Function is a two-dimensional function locally given by , , . The Jacobian matrix of is the constant matrix

To determine the matrix, one writes

where and . The first linear system gives and . Substituting these expressions in the second linear system we obtain

The superficial energy on triangle for the skin () is then given by

and the whole superficial energy is approximated by

For and we have


4.2 Numerical approximation of the discrete energy minimiser

For a given vector and the boundary conditions, we deduce vector , hence the continuous linear piecewise function we use to compute the discrete energy functional. We then build the operator

where is a given set of parameters. The numerical solution we seek provides the vector which minimises the energy of the discrete mechanical system. Conjugate gradients method is employed to determine the minimiser of the discrete functional .

4.3 Cost function discretisation

Let be a set of parameters. We deduce the discrete gravity free configuration which provides the operator . Then we compute the solution minimising the discrete energy functional . We deduce the final discrete breast after applying the deformation

with being the function that provides the new position. We assess the measures on and produce vector . In conclusion, we have define the discrete measures operator . On the other hand, the discrete cost function reads


We apply the iterative process to the discrete versions to get the best parameter set by calculating successive variation of the parameters () until for a given threshold .

5 Synthetic Cases

To assess the model quality and robustness, several numerical experiences are carried out using a manufactured solution. We solve the direct problem by prescribing a given set of parameters and compute the associated measurements and test the method ability to recover the initial parameters independently of the initial guess. To this end, we create a numerical breast setting the parameters and compute the three configurations in the gravitational fields (see Figure 4 for the case stand up) that provides the reference set of measurements

depending on the mesh characteristic size that corresponds to the mesh with 2948 vertices and about 12000 tetrahedrons.

Figure 4: Breast under gravity with the given set of defined parameters .

We want to evaluate the importance of the initial guess in the final results with initial guesses that went from 10% (approximate guess), 30% (reasonable guess), 60% (bad guess) and over from the reference parameters. The results (in mean) can be seen in table 1.

Parameter Loss
Granularity & Parameters Geometrical Mechanical
0,08% 0,39% 0,42% 0,89% 0,27% 2,77%
0,10% 0,47% 1,39% 0,79% 0,36% 2,40%
0,05% 1,20% 1,12% 0,79% 1,16% 2,32%
0,06% 0,98% 1,07% 0,76% 0,95% 2,28%
0,08% 0,68% 0,92% 0,65% 0,85% 2,14%
0,07% 0,72% 0,97% 0,62% 0,86% 2,13%
Table 1: Values of the parameter loss (in %) with different sized meshes after 25 iterations. For each mesh size, we selected 10 tests and the values shown are the means from those tests.

These results show the effectiveness of the algorithm in terms of converging towards the reference parameters with a value of loss close and in some cases inferior to with the exception of the parameter with a loss value around .

5.1 Robustness

We evaluate the computational effort of the method, particularly when increasing the number of vertices. We report the mean squared error (residual) between the reference measurements and the measurements obtained with different mesh sizes with the reference parameters (10 cases per mesh size). We also report the relative running time with respect to the coarser mesh, in function of the mesh size in table 2.

Mesh size 226 572 1014 1404 2127
Residual 1.52e-3 2.72e-4 7.26e-5 3.91e-5 9.13e-7
running time 1 3.4 7.8 11.2 20
Table 2: Measurement’s difference residual between the reference measurements and the measurements obtained with courser meshes with the same values for the breast parameters and the work effort which is spent by the optimisation process

Notice that the measurements obtained with the reference parameters and a mesh with just 572 nodes are approximately 5% different when comparing with the reference measurements. This means that by using coarser meshes we are, synthetically, introducing error to the measures.

We evaluate the precision for different size meshes (namely medium, medium-thin, thin and very-thin corresponding to a number of 572, 1014, 1404, and 2127 nodes respectively) and we report the values in figure 5.

Figure 5: Parameter evolution of the optimisation process for different meshes and comparison with the reference values. Top panel corresponds to the coarse mesh of 572 vertices while bottom panel is obtained with 2127 vertices. The intermediate panels correspond to the meshes with 1014 and 1404 vertices.

We observe the volatility for some parameters during the first estimation steps but we obtain stable approximations after iteration 15. Tables

3 and 4

give the mean and standard deviation value of each parameter for those tests in the interval between iteration 15 and 25.

Parameter Values (Mean)
Granularity & Parameters Geometrical Mechanical
0,35% 6,26% 48,70% 8,71% 2,17% 20,33%
0,24% 4,30% 30,18% 6,82% 2,19% 14,74%
0,09% 4,25% 17,93% 5,57% 1,64% 14,30%
0,05% 2,53% 4,21% 2,23% 1,20% 5,35%
Table 3: Parameter value difference (in %) between the reference parameters (

) and the parameter values obtained with courser meshes after the optimisation process with an initial guess with a variance up to

Parameter Values (Standard Deviation)
Granularity & Parameters Geometrical Mechanical
0,27% 1,83% 4,55% 7,22% 5,01% 21,15%
0,14% 0,64% 3,66% 4,49% 3,53% 12,19%
0,05% 0,31% 2,94% 2,52% 1,68% 7,10%
0,05% 0,23% 2,89% 1,53% 1,34% 4,01%
Table 4: Parameter value standard deviation (in %) between the reference parameters and the parameters estimated with courser meshes

Three levels/rates of convergence are highlighted:

  • High convergence. The geometrical parameters quickly converge to the solution after few iterations, and present a very stable behaviour (standard deviation values inferior to 2%).

  • Mild convergence. The estimation method presents good approximations for the physical parameters and .

  • Rough convergence. Parameters and approximation is far away when using coarse meshes and converge to a wrong solution after some few iterations.

5.2 Sensitivity Analysis

Sensitivity analysis is an important issue for assessing the robustness of the method and evaluating the measures that preponderantly influence the parameters. Moreover, potential errors might derive from the performed measurements in a practitioner’s office and one has to assess the consequences of such errors on the parameters. At last, it is a good indicator of ill-posed inverse problem. Minimising the cost function provides the relation , i.e. the best parameters that fit the measurements. The sensitivity analysis aims at assessing the impact of the measure variations on the parameters evaluation. Low sensitivity of parameter with respect to measure means that the relative partial derivative

while an high sensitivity is obtained if .

We assume that practitioners provide measurements with a relative error of . To simulate the error impact, we define two new sets of parameter and corresponding to a positive and a negative variation of of parameter , the others being fixed. To calculate the sensibility value we approximate the relative derivative with the expression


It results a matrix that approximate the Jacobi matrix. Computations are achieved with the very-thin mesh and we report the results in Table 6 with the parameters in row and the measures in columns.

Such data have to be handled and interpreted with caution since the variation derives from the we prescribe and computational errors due to the discretization may degrade the approximation accuracy. Therefore we have two variations to take into consideration

which corresponds to the error of the variation of the measure and

which corresponds to the error of the discretization with respect to the exact parameters

. We classify the data in three groups. The red groups are the data where the discretization error is preponderant, namely

. It results that the derivative approximation is wrong since the errors are mainly due to the discretization. The yellow corresponds to the data where , i.e. the discretization errors and the variations are of same order. At last, the green group provides the data such that which guarantee that the derivative is the consequence of the variation of parameter and can be considered as a correct approximation of the partial derivative. The values in the table are mapped according to the colours of reliability.

Figure 6: Table with the sensibility analysis of the parameters. (Red - unreliable values, yellow - slightly reliable values, green - reliable values)

Measures , and corresponding to the volume for the three positions have no significant impact on the evaluation of the parameters since the error deriving from the discretization are too large keeping from evaluating an approximation of the derivative. Such measurements should be discarded from the cost function. In the same way, the given by the measure is not relevant and should be eliminated. The other measures are valid since they provide an admissible approximation of the partial derivatives. We sum up hereafter the main conclusions.

  • Parameters and do not present a significant sensitivity to the measurements except with measurement corresponding to the in the lying down position.

  • We observe a great sensitivity for and with respect to the measurements that act in opposite way since we have positive and negative derivatives.

  • Parameter and have a high sensitivity with respect to the measurements over the 3 positions with positive or negative derivative leading to antagonist effects.

6 Conclusion

We have proposed a simple model to simulate breast displacement and identify the mechanical parameters associated to the Neo-Hookean model. The main points are the introduction of a two-parameter geometry for the gravity-free configuration and a user friendly set of measurements that does not require sophisticated equipment and turns to be adequate for routinely operations. The sensitivity study with respect to the measure and the mesh has been carried out in a synthetic context to demonstrate the robustness of the method and its capacity to retrieve the good parameters with a controlled range of error.

7 Conflict of Interest

The authors declare that there is no conflict of interest regarding the content of this article.


  • Affagard et al. (2015) Affagard, J. S., Feissel, P., and Bensamoun, S. F. Identification of hyperelastic properties of passive thigh muscle under compression with an inverse method from a displacement field measurement. Journal of Biomechanics, 48(15):4081–4086, 2015.
  • Azar et al. (2002) Azar, F., Metaxas, D., and Schnall, M. Methods for modeling and predicting mechanical deformations of the breast under external perturbations. Medical Image Analysis, pages 1–27, 2002.
  • Ball (1983) Ball, J. Energy-minimizing configurations in nonlinear elasticity. Proceedings of the International Cogress of Mathematicians, Warzawa, 1983.
  • Cardoso et al. (2003) Cardoso, A., Coelho, G., Zenha, H., Sá, V., Smirnov, G., and Costa, H. Computer simulation of breast reduction surgery. Aesth. Plast. Surg., pages 68–76, 2003.
  • Cardoso et al. (2010) Cardoso, A., Costa, H., Sá, V., and Smirnov, G. On the importance of chassaignac’s space in breast modelling. IV European Conference on Computational Mechanics, Paris, pages 16–21, 2010.
  • Ciarlet (1998) Ciarlet, P. Three-dimensional elasticity. Elsevier, Amsterdam, 1998.
  • del Palomar et al. (2008) del Palomar, A. P., Calvo, B., Herrero, J., López, J., and Doblaré, M. A finite element model to accurately predict real deformations of the breast. Medical Engineering & Physics, pages 1089–1097, 2008.
  • Eder et al. (2014) Eder, M., Raith, S., Jalali, J., Volf, A., and Settle, M. Comparison of different material models to simulate 3-d breast deformations using finite element analysis. Annals of Biomedical Engineering, 42(4):843–857, 2014.
  • Eiben et al. (2016) Eiben, B., Vavourakis, V., Hipwell, J. H., Kabus, S., Buelow, T., Lorentz, C., Mertzanidou, T., Reis, S., Williams, N. R., Keshtgar, M., and Hawkes, D. J. Symmetric biomechanically guided prone-to-supine breast image registration. Annals of Biomedical Engineering, 44(1):154–173, 2016.
  • Fung (2012) Fung, Y. Biomechanics - mechanical properties of living tissues. Springer-Verlag, 2nd edition, 2012.
  • Gascón et al. (2014) Gascón, B. H., Espés, N., na, E. P., Pascual, G., Bellón, J., and Calvo, B. Computational framework to model and design surgical meshes for hernia repair. Computer Methods in Biomechanics and Biomedical Engineering, pages 1071–1085, 2014.
  • Groves et al. (2013) Groves, R. B., Coulman, S. A., Birchall, J. C., and L.Evans, S. An anisotropic, hyperelastic model for skin: Experimental measurements, finite element modelling and identification of parameters for human and murine skin. Journal of the Mechanical Behavior of Biomedical Materials, 18:167–180, 2013. doi:
  • Han et al. (2012) Han, L., Hipwell, J. H., Tanner, C., Taylor, Z., Mertzanidou, T., Cardoso, J., Ourselin, S., and Hawkes, D. J. Development of patient-specific biomechanical models for predicting large breast deformation. Physics in Medicine and Biology, 2(57):455–472, 2012.
  • Isvilanonda et al. (2016) Isvilanonda, V., Iaquinto, J. M., Pai, S., Mackenzie-Helnwein, P., and Ledoux, W. Hyperelastic compressive mechanical properties of the subcalcaneal soft tissue: An inverse finite element analysis. Journal of Biomechanics, 49(7):1186–1191, 2016.
  • Lee and Zhuang (2012) Lee, H. and Zhuang, H. Biomechanical study on the edge shapes for penetrating keratoplasty. Computer Methods in Biomechanics and Biomedical Engineering, pages 1071–1079, 2012.
  • Lee et al. (2018) Lee, T., Vaca, E. E., Ledwon, J. K., Bae, H., Topczewska, J. M., Turin, S. Y., Kuhl, E., Gosain, A. K., and Tepole, A. B. Improving tissue expansion protocols through computational modeling. Journal of the Mechanical Behavior of Biomedical Materials, 82:224–234, 2018. doi:
  • Lorenzen et al. (2003) Lorenzen, J., Sinkus, R., Biesterfeldt, M., and Adam, G. Menstrual-cycle dependence of breast parenchyma elasticity: estimation with magnetic resonance elastography of breast tissue during the menstrual cycle. Investigative Radiology, 4(38):236–240, 2003.
  • Payan (2012) Payan, E. Y. Soft tissue biomechanical modeling for computer assisted surgery. Springer, Heidelberg, 2012.
  • Picinbono et al. (2003) Picinbono, G., Delingette, H., and Ayache, N. Non-linear anisotropic elasticity for real-time surgery simulation. Graphical Models, pages 305–321, 2003.
  • Poplack et al. (2004) Poplack, S., Paulsen, K., Hartov, A., Meaney, P., Pogue, B., Tosteson, T., Grove, M., and Soho, S. W. W. Electromagnetic breast imaging: Average tissue property values in women with negative clinical findings. Radiology, 2004.
  • Rajagopal et al. (2004) Rajagopal, V., Nielsen, P., and Nash, M. Development of a three-dimensional finite element model of breast mechanics. IEEE-EMBS 26th Ann. Intl. Conf., San Francisco, pages 5080–5083, 2004.
  • Rajagopal et al. (2007) Rajagopal, V., Chung, J.-H., Bullivant, D., Nielsen, M. F., and Nash, M. P. Determining the finite elasticity reference state from a loaded configuration. International Journal for Numerical Methods in Engineering, 72(12):1434–1451, 2007.
  • Zhong et al. (2005) Zhong, H., Wachowiak, M., and Peters, T. A real time finite element based tissue simulation method incorporating nonlinear elastic behavior. Computer Methods in Biomechanics and Biomedical Engineering, pages 177–189, 2005.