The solution of inverse non-linear elasticity problems that arise when locating breast tumours

Non-linear elasticity theory may be used to calculate the coordinates of a deformed body when the coordinates of the undeformed, stress-free body are known. In some situations, such as one of the steps in the location of tumours in a breast, the coordinates of the deformed body are known and the coordinates of the undeformed body are to be calculated, i.e. we require the solution of the inverse problem. Other than for situations where classical linear elasticity theory may be applied, the simple approach for solving the inverse problem of reversing the direction of gravity and modelling the deformed body as an undeformed body does not give the correct solution. In this study, we derive equations that may be used to solve inverse problems. The solution of these equations may be used for a wide range of inverse problems in non-linear elasticity.


Introduction
The accurate location of cancerous tumours in patients undergoing surgery is of considerable clinical importance.Tumours in women's breasts are often located using magnetic resonance (MR) imaging.When using this technique, the woman lies in the prone position with her breasts hanging downwards into the machine used to detect the signal.Imaging techniques are then used to locate the position of the tumour.Should the woman require surgery, then during the surgery she will usually lie in the supine position, and so her breast will occupy a different shape and position.The location of the tumour when the woman lies in this different position is not obvious, and so a technique for locating this tumour is, therefore, of much clinical interest.
The deformations of the breast described above are much too big for the theory of linear elasticity to be used, and so the theory of non-linear elasticity must be invoked.Non-linear elasticity theory relates the coordinates of an undeformed, stress-free body to the coordinates of the deformed, stressed body via a partial differential equation.When solving the forwards problem, the coordinates of the undeformed, stressfree body are known, and the governing partial differential equation may be used to calculate the coordinates of the deformed, stressed body.In particular, a numerical solution of the forwards problem may be computed using standard numerical techniques for solving partial differential equations.The procedure described above for locating the tumour and performing surgery requires that the breast takes two different deformed positions.In both cases, the deformation is independent of time, and caused by gravity.To map the location of the tumour from the prone position to the supine position requires two steps.The first step is to map the deformed breast in the prone position to a stress-free position, the position the breast would occupy in the hypothetical situation in which there is no gravity.The second step is to map the breast from the stress-free position to the position it will occupy when the woman lies in the supine position.The second of these steps is relatively easy if the constitutive relations for different tissue types are known, because it requires only the solution of a forwards problem as described above.However, the first step-the solution of the inverse elasticity problem-is more difficult.This is because for this step the coordinates of the deformed body are known, and we wish to calculate the coordinates of the undeformed body.We, therefore, know the solution of the governing partial differential equation for the corresponding forwards problem, but we do not know the computational domainthe region occupied by the undeformed body-on which this solution is defined.
In this study, we formulate equations that may be used to solve the inverse problem.We also demonstrate that, in general, the solution of inverse problems in non-linear elasticity is not as simple as reversing the direction of gravity and solving a forwards problem on the deformed body.Whilst the motivation of this work is the solution of inverse problems related to the location of tumours in breasts, the technique described is generic and may be applied to other problems in non-linear elasticity.One set of problems that could employ a technique for the solution of an inverse elasticity problem is image registration problems, where a deformable body organ is imaged in two different stressed positions.For example, X-ray images of the lungs are usually taken with the subject erect.However, the X-ray image may also be taken with the subject lying on their side should the clinician wish to investigate something that may show up better with the subject in that position.In addition, some subjects are so ill that it is only possible to obtain an X-ray image with the subject lying down.A soft tissue model of the lungs and the ability to solve inverse elasticity problems would aid in the fusion of data from these two images.A further set of examples is the determination of parameters for the stressstrain relationship of a given material.Gravity will always be acting on a material, and so it is impossible for a material to be in a truly stress-free state.The determination of parameters for the stress -strain relationship are often derived by loading the material with a known force.Stressstrain relationships may be evaluated by solving an inverse problem to estimate the stress-free configuration from the configuration due to gravity, followed by a forwards problem to calculate the configuration caused by loading the material.

Methods
Throughout this study, we use the summation convention, whenever, an index is repeated in a product, summation over all values of that index is implied.We define V 0 to be the region occupied by the undeformed body, and V to be the region occupied by the deformed body.These regions have boundaries ›V 0 and ›V, respectively.We begin by defining the stress and strain tensors that we use, and then write down the formulations of both the forwards and inverse problems in terms of these tensors.Finally, we demonstrate that the forwards and the inverse problem are identical in the special case, where quadratic and higher terms in the strain may be neglected, i.e. under the conditions where classical linear elasticity theory may be applied.

Stress and strain
Suppose X 1 , X 2 , X 3 are the coordinates of the undeformed body and x 1 , x 2 , x 3 are the coordinates of the deformed body.Then, the deformation gradient tensor F ¼ ðF iM Þ is defined by The Lagrange -Green strain tensor E ¼ ðE MN Þ is then defined in terms of the deformation gradient tensor by where d MN is the Kronecker delta [1].
The stress tensor can be defined in terms of the coordinates of the deformed or the undeformed body.We use two stress tensors in this study, the second Piola-Kirchoff stress tensor and the Cauchy stress tensor.The second Piola-Kirchoff stress tensor, T ¼ ðT MN Þ; is the force per unit undeformed area acting on the undeformed body.The Cauchy stress tensor, s ¼ ðs ij Þ is the force per unit deformed area acting on the deformed body.
For most materials a strain energy function, W, exists [2].Typical strain energy functions for incompressible materials are for biological tissue for Mooney-Rivlin materials where a and b are constants, p is the pressure (allowed to vary within the body) and the quantities I 1 , I 2 and I 3 are independent of any rotation of the coordinate axes and are given by [3] In this study, we shall restrict ourselves to simulations using the strain energy function for biological tissue, equation (3), although the technique which we describe is applicable to any valid strain energy function.The second Piola -Kirchoff stress tensor is then given in terms of the strain energy function by [3] The Cauchy stress tensor can then be calculated using the following relation between the Cauchy and second Piola -Kirchoff stress tensors [3] s

The forwards problem
When solving the forwards problem, we know the coordinates of the undeformed body, X 1 , X 2 , X 3 and wish to calculate the coordinates of the deformed body x 1 , x 2 , x 3 .We, therefore, formulate the forwards problem as a differential equation with X 1 , X 2 , X 3 as the independent variables and x 1 , x 2 , x 3 as the dependent variables.

The differential equation. The governing equations for static incompressible elasticity are
where r is the density of the body and g j is the gravitational force per unit mass in coordinate direction X j .We see, on using equations ( 2)-( 4) and ( 8), that equations ( 10) and ( 11) are differential equations with X 1 , X 2 , X 3 as the independent variables and x 1 , x 2 , x 3 as the dependent variables, as required.Note that equations ( 10) and ( 11) are four equations for four unknowns, three spatial coordinates x 1 , x 2 , x 3 and the pressure p.

Boundary conditions.
The boundary of the undeformed body is partitioned into two non-intersecting sets.On the first set, ›V D 0 ; there are Dirichlet displacement boundary conditions.On the second set, ›V ND 0 ; there are non-Dirichlet traction boundary conditions.The boundary conditions are T MN F jN N M ¼ s j j ¼ 1; 2; 3 on ›V ND 0 ð13Þ where ðx 0 1 ; x 0 2 ; x 0 3 Þ are the displacement boundary conditions, (N 1 , N 2 , N 3 ) is the outward pointing unit normal vector to the undeformed body and (s 1 , s 2 , s 3 ) is the force per unit undeformed area acting on the body.

The inverse problem
When solving the inverse problem the coordinates of the deformed body, x 1 , x 2 , x 3 are known and we wish to calculate the coordinates of the undeformed body X 1 , X 2 , X 3 .We, therefore, require a governing equation with x 1 , x 2 , x 3 as the independent variables and X 1 , X 2 , X 3 as the dependent variables.
2.3.1.The differential equation.We begin by writing the governing equations for static incompressible elasticity in terms of the Cauchy stress tensor [1]: However, equations ( 2), ( 3), ( 8) and ( 9) express s ij in terms of the quantities ›x i /›X M , and so we have not yet achieved our objective of writing a differential equation for X M ¼ X M ðx 1 ; x 2 ; x 3 Þ: Instead, we must first re-write the deformation gradient tensor, defined in equation ( 1), in terms of the quantities ›X M /›x i .Noting that we have and so we may write because, by equation (11), det F 21 ¼ det F ¼ 1: We now have F, and therefore (on using equations ( 2), ( 8) and ( 9)) E, T and s, in terms of the quantities ›X M /›x i .Equations ( 14) and ( 15) now become differential equations for the four dependent variables X 1 , X 2 , X 3 and p, as functions of the independent variables x 1 , x 2 , x 3 as required.We may now write equations ( 14) and ( 15) as Equations ( 19) and (20), with entries of F given by equation (18), are therefore, the equations that govern the solution of the inverse problem.
2.3.2.Boundary conditions.We partition the boundary of the deformed body into two sets in the same way as for the undeformed body.On the first set, ›V D , there are Dirichlet displacement boundary conditions.On the second set, ›V ND , there are non-Dirichlet traction boundary conditions.The boundary conditions are where ðX 0 1 ; X 0 2 ; X 0 3 Þ are the displacement boundary conditions, (n 1 , n 2 , n 3 ) is the outward pointing unit normal vector to the deformed body and (t 1 , t 2 , t 3 ) is the force per unit deformed area acting on the body.

The necessity of equations (19) and (20)
In classical linear elasticity theory, very few workers distinguish between the second Piola -Kirchoff stress tensor and the Cauchy stress tensor.As a result, properties of the original configuration of the body, such as, whether it is in the stressed or stress-free state, are rarely taken into consideration.If the stressed body is treated as a stressfree body, the formulation of the inverse problem is a forwards problem where the coordinates of the deformed and undeformed body have been interchanged and the direction of gravity has been reversed.In this section we investigate the validity of using this approach.
Suppose we were to attempt to solve the problem in the manner described in the previous paragraph.Then we may solve equations ( 10) and (11), interchanging the x i 's and X M 's, interchanging V and V 0 , and reversing the direction of gravity.Equation ( 11) is unchanged, equation (10) becomes which may be written where and The strains used in this calculation would be those calculated using equation (25), and so the strain energy function Ŵ for this problem is defined by On using equation ( 16), we may write and so equation ( 23) may be written which is different to the true governing equation of the inverse problem, equation (19).Hence, in general, this approach to solving the inverse problem is not valid.
2.4.1.Approximate solution when quadratic and higher order terms in the strain may be neglected.Writing the displacements as x i ¼ X i þ u i ; i ¼ 1; 2; 3 we may write the deformation tensor as Suppose the displacements are so small that we may neglect quadratic and higher order terms in u i , as is the case in classical linear elasticity.Define When neglecting quadratic and higher order terms in u i , we see from equations ( 30) and (31) that we may write the inverse of F as and we may then use equations ( 24) -(28) to give we see that, on neglecting quadratic and higher order terms in the strain, equation ( 19) is identical to equation (29).We, therefore, see that under conditions where quadratic and higher order terms in the strain may be neglected, i.e.where the assumptions in linear elasticity theory are valid, there is no need to distinguish between the stress-free and the stressed body.

Numerical methods for soft tissue modelling
Previous work on soft tissue modelling has been applied to the breast [5][6][7][8][9], the heart [10][11][12] and the liver [18].All of these authors use the finite element method to solve the governing equations although, as full details of the scheme used are sometimes not given, it is not always possible to comment on the suitability of the scheme used.
Breast tissue is assumed to be incompressible by most authors [5 -9].However, with the exception of [8], other workers do not enforce incompressibility using the mathematically correct procedure described above, where the constraint equation (11) and the additional dependent variable pressure are added to the governing equation for non-linear elasticity, equation (10).Instead, incompressibility is modelled by setting the Poisson ratio of breast tissue to 0.5 2 e, where e is a small, positive number.Having enforced incompressibility in this manner, a variety of methods were used to solve the governing equations.Schnabel et al. [9] used linear elasticity to approximate the displacement of breast tissue.Azar et al. [5,13] model the whole displacement as consisting of a large number of tiny displacements that are computed using linear elasticity, this approach allows an exponential stress -strain relationship to be built into the model.Samani et al. [6] and Liu et al. [7] used non-linear elasticity when computing displacements, allowing more general stress -strain relationships to be built into the model.
There are numerical complications caused by modelling incompressible tissue as the limit of the Poisson ratio tending upwards to 0.5.Standard finite element schemes, where the displacement is approximated by a continuous piecewise linear function, exhibit the phenomenon known as locking [14,15].Locking is a problem caused by these schemes being incapable of resolving both incompressibility and traction conditions, and leads to poor convergence rates as the mesh is refined.Locking may be avoided either by using schemes where the continuity requirement is relaxed, or schemes with a higher order polynomial approximation to the displacement.In all the publications on breast modelling cited above where authors described the finite element scheme used, the scheme used continuous piecewise linear functions for the approximation to the displacement.Nash and Hunter [12], who model passive cardiac tissue as being incompressible, avoid locking by (i) using a Hermite cubic approximation to the displacements, and (ii) modelling compressibility as has been described in section 2.2.

Numerical simulation
In this section we demonstrate numerically that equation (23) does not give the true solution to inverse problems in non-linear elasticity, and that equations ( 19) and (20) should instead be used to calculate the true solution.

Description of simulation
For ease of presentation we simulate a displacement in two dimensions, although qualitatively similar behaviour is seen in three dimensions.We simulate the displacement of the unit square.The edge y ¼ 0 is fixed, all other boundaries are free to move and are stress-free.We use the strain energy function for biological tissue given in equation ( 3) with parameters a ¼ b ¼ 1; this strain energy function is used for isotropic tissue with an exponential relationship between stress and strain.Density was chosen so that r ¼ 1: Gravity was given by g ¼ ð0; 2Þ t : We fully accept that this simulation is not representative of deformations of a breast, this simulation is in only two dimensions; the shape is different; no account is taken of different tissue types; the material is assumed to be isotropic; and the parameters used have no physiological significance.In reality, a breast will consist of at least three tissue types, fibroglandular tissue, fatty tissue and skin.These different tissues will all have different stress -strain relationships, different degrees of anisotropy and different densities.However, the derivation of equations ( 19) and (20) permitted the use of a very general strain energy function.This algorithm, therefore, allows us to use many tissue types, anisotropic tissue and a general computational domain.We have chosen this simulation purely on the basis that a numerical solution may be calculated using standard numerical techniques.This simulation is included in this study only to demonstrate that, in general, using equation ( 23) is not a suitable technique for solving inverse problems in non-linear elasticity and that equations ( 19) and (20) should be used instead.
We denote the coordinates of the unit square by X ¼ ðX 1 ; X 2 Þ: For the simulation described above, we perform three calculations.
(1) We calculate the coordinates of the deformed body, x ¼ ðx 1 ; x 2 Þ; as a function of X using equations ( 10) and ( 11).( 2) Using the values of x from the first calculation, we solve the equations that govern the inverse problem, equations ( 19) and ( 20). ( 3) Using the values of x from the first calculation, we solve equations ( 23) and (11) to calculate the predicted coordinates of the undeformed body by reversing gravity and making no distinction between the stressed and stress-free body, as described in section 2.4.We denote these predicted coordinates by X: The maximum displacement of the forwards problem, D 0 , is given by We then define D 1 (X) to be the distance between the true deformed coordinates and the undeformed coordinates predicted by equations ( 23) and (11): As an indicator of the error induced by using the incorrect formulation of the inverse problem we calculate the maximum and the average value of D 1 , scaled by the maximum displacement of the forwards problem, D 0 .We therefore calculate the following quantities: where D 1 is the area weighted average of D 1 , given by

The finite element solution
The governing equations were solved using the finite element method, see for example, [16].The unit square was discretised into 800 uniformly sized elements, there were 20 elements in the x direction and 40 elements in the y direction.A bilinear approximation was used for the dependent spatial variables in each of these elements.For a stable approximation to the pressure as element size is reduced, a lower order approximation must be used to calculate the pressure, see for example, [16].We, therefore, use a piecewise constant approximation on each element to calculate pressure.The non-linear equations arising in the finite deformation calculations were solved using Newton's method (see for example, [17]).

Results of the simulation
In figure 1(a) we plot, the boundary of the deformed body (solid line); the boundary of the region that is the solution of the inverse problem (dashed line); and the boundary of the region that is calculated by reversing the direction of gravity (dot-dashed line).In figure 1(b), we show a close up of figure 1(a), to show more detail of the region around x ¼ 0: It is clear from figure 1 that the true solution of the inverse problem (dashed line) differs significantly from the solution of the inverse problem that has been calculated by reversing the direction of gravity (dotdashed line).This is confirmed in table 1 where we list the errors that are induced by the incorrect formulation of the inverse problem.We see that the maximum error is 23% of the true displacement.The magnitude of the error seen here would not be acceptable for the location of a tumour.

Discussion
We have seen in this paper that, unless quadratic and higher order terms in the strains may be neglected, the solution of inverse problems in non-linear elasticity is more complex than simply reversing the direction of gravity, interchanging the role of the deformed and undeformed coordinates, and solving the resultant forwards problem.We have therefore formulated equations (equations ( 19) and ( 20)) that express the coordinates of the undeformed body in terms of coordinates of the deformed body, and may therefore be used to solve the inverse problem.Solving these equations allows us to solve the inverse elasticity problem described in the Introduction, which was one of the steps that are used to locate tumours in breasts.In addition, the technique is valid for any suitable strain energy function, and so may be applied in other areas where the solution of inverse elasticity equations is required.
It should be remembered that many assumptions are made when deriving a mathematical model of soft tissue.The assumptions made for the simulations in this study are listed in section 3.1.These assumptions will all induce errors in the computed displacements to some degree.For example, the fact that motivates this study-the location of tumours in breasts-each individual breast is of different shape, consists of different tissue types, and has a different stress-strain relationship.Quantifying the error induced by making these assumptions is very hard to predict.However, it is evident that in situations where non-linear elasticity is an appropriate mathematical model, the solution of inverse problems requires careful consideration.In these situations solving inverse problems correctly eliminates one component of the total error-the error that would be induced by solving a forwards problem with the direction of gravity reversed instead of solving the inverse problem.
The author would like to thank Professor Peter Hunter, Dr Martyn Nash and Dr Poul Nielsen of the University of Auckland for hospitality and many enlightening discussions during the preparation of this manuscript.We also thank Professor Sir Mike Brady and Professor David Gavaghan of the University of Oxford.The author is pleased to acknowledge the financial assistance of the

Figure 1 .
Figure 1.The boundary of the deformed body (solid line), the boundary of the region calculated by solving the inverse problem (dashed line) and the boundary of the region that is calculated by reversing the direction of gravity (dot-dashed line).The whole simulation is shown in (a) and (b) shows a close up of the region around x ¼ 0: See text for further details of the simulation.

Table 1 .
Calculated values of D 0 , max(D 1 )/D 0 and D 1 =D 0 for the simulation described in the text.