Elastography Method for Reconstruction of Nonlinear Breast Tissue Properties

Elastography is developed as a quantitative approach to imaging linear elastic properties of tissues to detect suspicious tumors. In this paper a nonlinear elastography method is introduced for reconstruction of complex breast tissue properties. The elastic parameters are estimated by optimally minimizing the difference between the computed forces and experimental measures. A nonlinear adjoint method is derived to calculate the gradient of the objective function, which significantly enhances the numerical efficiency and stability. Simulations are conducted on a three-dimensional heterogeneous breast phantom extracting from real imaging including fatty tissue, glandular tissue, and tumors. An exponential-form of nonlinear material model is applied. The effect of noise is taken into account. Results demonstrate that the proposed nonlinear method opens the door toward nonlinear elastography and provides guidelines for future development and clinical application in breast cancer study.


Introduction
Breast cancer is one of the major threats to public health all over the world. Currently, X-ray mammography is the primary method for early detection and characterization of breast tumors [1]. While more efficient in detecting malignancies as age increases or the breast becomes fatty, mammography fails to detect small cancers in dense breasts. Further, mammography may not be specific in terms of tumor benignity and malignancy [2][3][4].
To solve these problems associated with mammography, a number of technologies have been explored. Detection and characterization of breast tumors can be enhanced by recognizing the difference of elastic modulus (stiffness) among normal soft tissues and malignant and benign tumors. Elastic properties of breast tissues may become an indicator of histological diagnosis [5]. An imaging technology called elastography was developed as an approach to imaging tissue elastic modulus in a quantitative manner for detection of breast tumors in 1990s. The general basis of elastography is to induce motion within tissues under investigation by either an external or internal mechanical stimulation. Conventional medical imaging modalities are then used to measure the spatial deformation, from which the mechanical properties can subsequently be reconstructed. Most simulations are either ultrasound or magnetic resonance (MR) imagebased [6][7][8][9][10] that take the dynamic or quasistatic interior displacement field, completely or partially, as input for identification of the elasticity properties. Current elastography reconstruction framework is based on the assumption of linear elasticity theory. It is shown, however, that the deformation of most biological soft tissues are not linearly elastic [5,11]. Consideration of nonlinear models is essential for elastography in clinical applications.
The present work aims at developing a nonlinear elastography model for breast tissues. We introduce a nonlinear adjoint gradient method that significantly improves the numerical efficiency and enhances the stability of elastography reconstructions. In Section 2, we present an algorithm based on nonlinear adjoint gradient method. We apply the iterative Newton method to solve unknown displacements and forces. In order to find elastic modulus distribution that 2 International Journal of Biomedical Imaging minimizes the objective function based on measured and calculated forces, the adjoint gradient method is employed to provide user-supplied gradients in nonlinear elastography. Numerical simulation is described in Section 3 including establishing a three-dimensional (3D) heterogeneous FEM breast phantom and applying exponential-form nonlinear material model. Four types of compressive loadings are applied in the forward problem. The iterative reconstruction based on force measurements on surface is also detailed. In Section 4, the result of the inverse problem is obtained and the effect of noise is investigated. A user-defined penalty function is introduced to reduce the impact of noise on reconstruction.

Finite-Strain Deformation Equations.
We use continuum description for the breast tissues. Let Ω 0 be a biological object. From the displacement tensor u(X) based on Lagragian coordinate system X, the deformation gradient is F = I + ∂u/∂X and the Green strain is E = (F T · F − I)/2 where "·" denotes the contraction operation between two tensors. The breast tissues are assumed hyperelastic so that the second Piola-Kirchhoff stress is S = ∂W(E; p)/∂E, where W is the strain energy and p denotes material parameters in the model. We use ";" to separate material parameters from deformation variables. The governing equation and boundary conditions for u are ∇ · S · F T + b(X, u(X)) = 0 X ∈ Ω 0 , N(X) · S · F T = t(X, u(X)) X ∈ Γ 0 t , The boundary of Ω 0 consists of Γ 0 t , with N the outer normal, on which external force t is applied, and Γ 0 u where displacement u is prescribed. Here, we consider general problems that the body force b and surface force t are deformation dependent. Following the standard finiteelement (FE) method [12], the displacement u is discretized as nodal displacement vector {u} = {u 1 , u 2 } T , where u 2 corresponds to u prescribed on Γ 0 u , and u 1 is to be solved from nonlinear equations: The internal nodal force f in corresponds to stress S; that is, it changes with u 1 and material parameters p, as u 2 is given. The external nodal force f out 1 is due to the prescribed surface force t on Γ 0 t and body force b in Ω 0 . It varies with the displacement in large deformation. The nodal force f out 2 is the unknown constraint force on Γ 0 u . A classic Quasi-Newton method is employed to solve (2) for u 1 as shown in Appendix A.

Nonlinear Elastography Algorithm.
Experimental measurements for elastography include displacement and force.
We consider that the biological object Ω 0 is discretized into FE mesh, and the measurements are discretized consistently into nodal displacement and nodal force. We catalog the measurements as the following. (i) If the displacement at a node is known (prescribed or measured after deformation), it will be included into u 2 which is considered "prescribed" in FE (2). The corresponding nodal force (known or unknown) will be in the constraint force, denoted as F M 2 . Note that F M 2 corresponds to f out 2 in (2). (ii) All the other nodal displacements will be in u 1 and the corresponding nodal force will be in f out 1 . For category (ii), f out 1 must be considered "prescribed" to fulfill the requirement of the wellpostness of a solid-mechanics problem.
For elastography problem, the obtained u 2 and f out 1 are known in (2), and the constraint force F M 2 is considered as "measurement." For given material parameters p, the unknown displacement u 1 and constraint force f out 2 (which depends on p) will be solved from FE (2). Elastography procedure thus searches for p so that the overall difference between measured F M 2 and computed f out 2 is minimum, that is, minimize objective function where the diagonal weight matrix Λ = diag(a 1 , a 2 , . . . , a j , . . .). The component a j = 1 when the jth component of F M 2 is measured or a j = 0 otherwise. The present algorithm is mathematically equivalent to our previous work [10] where the displacement is used as "measurement." For breast tissue whose tangent stiffness significantly increases with strain, this "force version" shows better numerical efficiency.
Efficient and robust optimization-based elastography schemes request user-supplied gradient ∂Φ/∂p. The previous elastography studies used direct or approximate finitedifference method [13,14] for the gradient calculation. The computational expense of these methods increases proportionally with the number of material parameters, and becomes unaffordable for problems involving finite-strain deformation. Recently, an adjoint method was introduced to compute the gradient analytically [10,[15][16][17]. The corresponding nonlinear finite element formulas are shown in Appendix B. After u 1 is solved for given p, the gradient is readily obtained as where the virtual adjoint displacements w 1 and w 2 are solved from linear equations with the tangent stiffness matrix K eff defined in (A.2). The most significant features of the adjoint method are its analytical formulation, high accuracy, and computational   efficiency [18]. Since K eff and its LU factorization have been calculated when solving the FE (2), as in (A.1), the additional expense for w 1 and w 2 in (5) is minimal. Furthermore, it only needs to solve one linear (5) regardless of the number of unknown parameters in p. As shown in Figure 1, a well-tested L-BFGS subroutine is applied for the present optimization-based nonlinear elastography.

Numerical Simulations
This section presents phantom simulations to identify the nonlinear elastic properties of the fatty, glandular, and cancerous tissues in a breast. First of all, a 3D breast FEM phantom attracting from the real data is introduced. Fung's model [19] is applied to describe the deformation of breast. With the applied loading, the forward-problem computation is performed. Furthermore, boundary forces are extracted from the forward computation results and are used as input for reconstruction for a breast phantom.

Breast Phantom.
To perform numerical simulations, a 3D numerical heterogeneous breast phantom extracting from real CT images, containing fatty tissue, glandular tissue, and a tumor is established (Figure 2). Boundaries of these regions are described with sets of splines. The phantom is discretized with standard 3D tetrahedral elements, consisting of 7303 elements and 1583 nodes ( Figure 3).
It is well known that the mechanical behavior of biological soft tissue is nonlinear. Hyperelastic models have commonly been used to represent the stress-strain relation of biological soft tissue [20,21]. Fung and coworkers developed a set of hyperelastic models for bio-tissues [20]. In this work, we employed a Fung-type isotropic strain-energy function for the breast tissues, as W(E; p) = γ(exp(λ(E : I) 2 )/2 + μE : E) − 1 , where E is the Green strain, p = {γ, λ, μ} denotes the material parameters. The second Piola-Kirchhoff stress S is thus Specifically, in uniaxial tension/compression, the relation between the axial strain E and axial stress S is where E Young = μ(3λ + 2μ)(λ + μ) −1 is the Young's modulus. Equation (7) is the exponentially nonlinear model that can be applied for breast tissues. Two parameters γ and E Young can be determined by the experiment [22]. Wellman [5] developed a technique to measure the nonlinear elastic parameters of breast tissues using force-displacement data of thin slice tissues undergoing indentation experiment. The tissue samples tested were obtained during surgery and were tested immediately after removal from the body. Six breast tissues (fatty, gland, lobular carcinoma, fibroadenoma, infiltrating ductal carcinoma, and ductal carcinoma in situ) were tested with their nonlinear stress-train curves shown in Figure 4. It is shown that the mechanical properties between normal soft breast tissue and tumors are quite different. Tumors are stiffer than the surrounding breast tissues and malignant tumors are much stiffer than benign ones.    (7)).

Forward Problem.
After the breast phantom and material model are established, a forward problem is solved in which material parameters and external loading are given and deformation is solved. Displacements u = 0 are prescribed on the base of phantom. Tilted compression by paddles is applied on upper surface ( Figure 5). The paddle close to tumor gives compression and another paddle is fixed during loading. Four types of compression loading with different angles are applied. As shown in Figure 6, blue lines represent paddle locations before loadings, green lines represent those after loadings. Figure 7 shows the comparison of paddle locations in four loadings. Note that the right paddle is fixed for all loadings.
The reason that four sets of loadings are applied in forward problem is to provide more information to reconstruct material parameters. Most of inverse problem in elasticity is nonuniqueness. Previous research [10] has demonstrated that one set of measurements of displacements and forces may not provide sufficient information of the reconstruction of modulus distribution.
Given material parameters and loadings, the displacements and forces can be calculated. The surface force will be used as input to reconstruction material parameters in the inverse problem. In fact, surface displacement and force are equivalent as input to solve inverse problem. Most of previously research applied displacement in linear elastography. In this nonlinear elastography study, however, it is found the reconstruction is more sensitive to force than displacement. Therefore surface force is measured and compared with calculated force to reconstruct material parameters.

Inverse Problem.
Reconstruction for nonlinear elastic moduli in 3D breast phantom take input extracted from the deformation in response to loading modes A∼D described in Figure 6. In each loading, the forces on surface nodes are measured as input. Then initial guesses for material parameters are given. In this study, the same initial guesses are applied to three materials: λ 0 = 20, μ 0 = 10, and γ 0 = 1. Therefore, the surface force can be calculated based on initial guesses of material parameters. The error between calculated and measured force is used to set up the objective function (3). Following the iterative optimization procedures (Figure 1), gradients of the objective function are calculated and material parameters are updated to approach the optimal values. International Journal of Biomedical Imaging

Ideal Input.
For three materials, elastic parameters are: λ f = 35, μ f = 12.5, and γ f = 0.4 for fatty tissue; λ g = 50, μ g = 25, and γ g = 0.25 for glandular tissue, and λ d = 80, μ d = 35, and γ d = 1.5 for ductal carcinoma in situ. Table 1 shows the initial estimate and reconstructed results for nonlinear elastography. The results in the first part are based on the ideal input. It is demonstrated that the reconstructed results are very close to the real values. The maximum error is 1.916% (λ for tumor) since the effect of the tumor on surface force measurement is the smallest. A similar situation occurs in linear elastography reconstruction.

Adjoint Methods.
Elastography includes forward and inverse problem. In forward problem, material parameters and loadings are given to calculate the deformation; while in inverse problem, external loadings and deformation are known to reconstruct material parameters. Most researchers established certain objective function and minimized it with a proposed iterative algorithm. The challenge is how to calculate the gradient of objective function efficiently and accurately. A straightforward calculation of gradients requires solving stiffness matrix in each iteration, which takes most of the time consumed in the finite element method. Table 1: Initial estimate and reconstructed results for nonlinear elastogrpahy. The results in the first part are based on ideal input (without noise). The ones in the second part are based on input with 5% noise and regularization is not used. The third part is based on input with 5% noise and regularization is applied to reduce the impact of noise. (λ and μ are dimensionless, γ is in kPa.)

Fatty
Glandular In this study an adjoint method is employed to analytically calculate the gradients. Oberai et al. [16] adopted an adjoint method and proposed a numerical scheme for reconstructing the nonuniform shear modulus field for incompressible isotropic materials using one component of displacement field. Liu et al. [10] applied this method for anisotropic materials. In this study, the adjoint method is applied for nonlinear elastography. The advantage of adjoint method is to solve two adjoint displacements during each iteration (w 1 and w 2 at (B.2) in Appendix B), instead of the whole stiffness matrix, that increases the numerical efficiency significantly. In [18] Oberai et al. compared three different iterative methods: (1) a gradient-based method where the adjoint approach is used to calculate the gradient; (2) a gradient-based method where the straightforward approach is used to calculate the gradient, and (3) the Gauss-Newton method. The results show that "leading-order costs for the gradient-based method with the adjoint approach are smaller than the other two methods." In fact, without the adjoint method, nonlinear elastography can only be discussed in concepts [23] or applied on simple objects using supercomputing power [24]. While in this study, material parameters in a real breast FEM phantom are reconstructed accurately by adopting this adjoint method. The results are encouraging for further clinical applications.

Multiple Sets of Measurements.
Above results are based on four sets of force measurements on surface. For 2D isotropic elastography, Barbone and Bamber [25] has shown that one set of measurements for the displacements and forces, especially those taken only on the boundaries, may not provide sufficient information to the reconstruction of modulus distribution, because of the nonuniqueness nature of most inverse problems in elasticity. Barbone and Gokhale [26] further demonstrated the feasibility of using multiple displacement fields to reduce the likelihood of nonuniqueness for 2D isotropic elastography. Liu et al. [10] discussed the multiple sets of measurements in 3D anisotropic media. However research effort for 3D nonlinear elastography has not been found in literatures. In this study four independent tilted compression loadings are designed and surface force is measured to reconstruct internal material parameters. Due to different initial guesses, stable material parameter can be reconstructed, showing that the multiple sets of measurements are feasible to obtain the unique solution of inverse problem in nonlinear elastography.
The key points for using multiple sets of measurements are to bring more deformation modes simultaneously into consideration. The loadings should be close to tumors in order to make tumors have larger deformation. In Figures 5-7, four sets of loadings close to tumors are designed to obtain more deformation for reconstruction. This explanation serves as a guideline for design of loadings in clinical applications. Reduction of the number of required loadings will increase the clinical efficiency and benefit the patients. In other words, loadings with the richest stressstrain information are most desired. The design of feasible loadings is important for success in clinical application of nonlinear elastography.

Input with Noise.
The above results are based on ideal input. However, noise cannot be avoidable in experiments. Its impact on reconstruction is therefore investigated with 5% of noise applied on surface measurements. The results are shown in the second part of Table 1. It is obvious that the reconstructed parameters are away from the real values. A possible reason is that, while the goal is to find real material parameters that minimize the objective function, a global minimum is not well defined with noise. In another words, real material parameters may not give a global minimum of the objective function with noise. Several similar parameters around real ones may give some local minimums. Our algorithm may reach one of them but fails to find real one because they all give local minimums. A typical way to solve this problem is to add a penalty term which provides additional constraint on the solution space. The penalty term may push the solution into the right area of the solution space and minimize the resulting objective. A new objective function is therefore proposed as where Φ is the original objective function (3) which represents the difference between measured and calculated values.
International Journal of Biomedical Imaging 7 In addition, χ is the regularization factor and Π is the penalty term. Specific forms of penalty term have been designed for different problems [27]. In this study, an exponential form of penalty term is applied as where K = 9 is the total number of material parameters, and ξ k and a k are the reconstructed and true elastic parameters, respectively. If the true material parameters are unknown, a k can be estimated as close as possible. This is reasonable since several experiments have provided the nonlinear elastic parameter for breast tissues.
The results are shown in the third part of Table 1, demonstrating that the regularization method significantly improves the reconstruction data. Most of parameters are close enough to the true values. The largest error occurs for χ in fatty tissue with approximately 15%. By adding the penalty term, the reconstructed values are pulled from local minimum into the range close to true values. It is noted that the regularization factor χ varies in different scenarios. For example, χ can be smaller for higher confidence of experiment-data accuracy. On the contrary, if a k is closer to the true value, χ should be larger. For this study, the elastic parameters in fatty and glandular tissues are stable, comparing with the ones in tumor (Figure 4). Therefore, larger regularization factors are used for fatty and glandular tissues while smaller χ can be applied for tumors.

Conclusions
This paper presents a study on nonlinear elastography of biomedical imaging, in which a 3D model is developed for heterogeneous breast tissues extracting from real images including fatty tissue, glandular tissue, and tumors. Based on the large-deformation constitutive law, discretized nonlinear equations are solved for displacement, strain, and stress fields in breast tissues with given tumors under external compression at breast boundaries. A 3D inverse-problem algorithm is developed to reconstruct the material parameters for nonlinear elastic constitutive relation of breast phantoms with tumors. The adjoint gradient method is introduced to improve the numerical efficiency and enhance the stability of elastography reconstruction. Results demonstrate that this work opens the door toward nonlinear elastography, and provides guidelines for future developments and clinical applications in breast cancer study.

A.
A classic quasi-Newton iterative method [28] is employed to solve (2) for u 1 . Let u 1 (n) be the trial solution of the unknown u 1 at the nth iterative step. An improved solution u 1 (n+1) = u 1 (n) + δu 1 can be obtained at the next step, in which δu 1 is the solution of linear equations: The iteration is put into effect until the residue f out 1 − f in 1 is smaller than the preset criterion.

B.
It is difficult to calculate gradient ∂Φ/∂p because u 1 is an implicated function of P. Here we will derive an adjoint method for the gradient ∂Φ/∂p in (4). We introduce the constraint (2) into the objective function (3), and obtain a Lagrangian: where w 1 and w 2 are arbitrary virtual displacements. It is noted that variables u 1 and p are no longer coupled. It is also noted that Φ = L and δΦ = δL under the constraint (2). Since the equality constraint (2) is satisfied with u 1 , the variation δL can be expressed as give (5) in the text. By introducing the adjoint method, it seems that more equations and variables (w 1 and w 2 ) need to be solved. But the solution of (B.6) is straightforward and computational cost is minimized because K eff = ∂( f in 1 − f out 1 )/∂u 1 has been computed and factorized when solving for the displacement u 1 as in Appendix A. For nonlinear problem, it is very expensive to calculate ∂u 1 /∂P, while the adjoint method avoid to solved ∂u 1 /∂P by introducing variables, w 1 and w 2 , which could be solve in simple linear equations.