1 Introduction
Biomechanical breast modelling and simulation with the finite element method (del Palomar et al., 2008, Eder et al., 2014) provides powerful analysis tools to predict 3D breast deformations for various biomedical applications. Breast tissue modelling is usually based on nonlinear mechanical models such as the Neohookean model expressed in the stressfree unloaded reference state of the breast.
Both the parameters of the mechanical model and the initial gravityfree 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 lowcost 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 fourparameters nonlinear mechanical model defining on a simplified twoparameters 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 stressfree 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 twoparameter 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 welldeveloped in the biomechanics 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 Neohookean is considered as a welladapted 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 stressfree geometrical domain of the breast equipped with the Neohookean 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 gravityfree 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 nontruncated 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 .
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 inframammary fold.
The hyperelastic neoHookean compressible relations (Cardoso et al., 2010) equipped with the socalled Chassaignac space represent a wellaccepted 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 neoHookean 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(1) 
subject to the constraints
(2)  
(3) 
where is the Chassaignac coefficient, and stands for the bulk density. Constraint (2) represents the fixation of the inframammary 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 strainenergy and the skin strainenergy 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 twoparameters 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 stressfree configuration , as shown in Figure 1, which defines the operator
(4) 
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
(5) 
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.
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
(6) 
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 GaussNewton method to take advantage that function casts in the specific form
(7) 
The firstorder Taylor series expansion reads
(8) 
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
(9) 
We seek a perturbation that minimises , i.e., . Therefore Equation 9 provides the expression
(10) 
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 nonoverlapping 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 subvector 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 twodimensional 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
and
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
(11) 
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.
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  
R  H  

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% 
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.52e3  2.72e4  7.26e5  3.91e5  9.13e7 
running time  1  3.4  7.8  11.2  20 
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, mediumthin, thin and verythin corresponding to a number of 572, 1014, 1404, and 2127 nodes respectively) and we report the values in figure 5.
We observe the volatility for some parameters during the first estimation steps but we obtain stable approximations after iteration 15. Tables
3 and 4give 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  
R  H  

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% 
) and the parameter values obtained with courser meshes after the optimisation process with an initial guess with a variance up to
fromParameter Values (Standard Deviation)  
Granularity & Parameters  Geometrical  Mechanical  
R  H  

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% 
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 illposed 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
(12) 
It results a matrix that approximate the Jacobi matrix. Computations are achieved with the verythin 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.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 NeoHookean model. The main points are the introduction of a twoparameter geometry for the gravityfree 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.
References
 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. Energyminimizing 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. Threedimensional 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 3d 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 pronetosupine breast image registration. Annals of Biomedical Engineering, 44(1):154–173, 2016.
 Fung (2012) Fung, Y. Biomechanics  mechanical properties of living tissues. SpringerVerlag, 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: https://doi.org/10.1016/j.jmbbm.2012.10.021.
 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 patientspecific 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., MackenzieHelnwein, 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: https://doi.org/10.1016/j.jmbbm.2018.03.034.
 Lorenzen et al. (2003) Lorenzen, J., Sinkus, R., Biesterfeldt, M., and Adam, G. Menstrualcycle 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. Nonlinear anisotropic elasticity for realtime 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 threedimensional finite element model of breast mechanics. IEEEEMBS 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.