Abstract
Quantification of the mechanical behavior of normal and cancerous tissues has important implication in the diagnosis of breast tumor. The present work extends the authors' nonlinear elastography framework to incorporate the conventional X-ray mammography, where the projection of displacement information is acquired instead of full three-dimensional (3D) vector. The elastic parameters of normal and cancerous breast tissues are identified by minimizing the difference between the measurement and the corresponding computational prediction. An adjoint method is derived to calculate the gradient of the objective function. Simulations are conducted on a 3D breast phantom consisting of the fatty tissue, glandular tissue, and cancerous tumor, whose mechanical responses are hyperelastic in nature. The material parameters are identified with consideration of measurement error. The results demonstrate that the projective displacements acquired in X-ray mammography provide sufficient constitutive information of the tumor and prove the usability and robustness of the proposed method and algorithm.
1. Introduction
Breast cancer is a major threat to public health in the world. In USA and Europe, approximately 10% of women develop breast cancer during the course of their lives. While the specific causes of breast cancer are unknown, early detection and characterization of breast tumors is the key to successful treatment. Currently, X-ray mammography, a low-dose X-ray imaging modality, is the primary diagnosis method in clinics [1]. While being more efficient in detecting malignancies as age increases or the breast becomes fatty, mammography fails to identify small cancers in dense breasts. Furthermore, mammography may not be specific in terms of tumor benignity and malignancy. About 80% of suspicious masses referred by mammography for surgical breast biopsy are in fact not malignant [2–4]. These false-positive mammograms may induce patients’ anxiety, distress, and intrusive thoughts.
A number of techniques have been attempted to address these problems associated with mammography. From the viewpoint of mechanics, the tissue stiffness is an important index for diagnosis of breast cancers, as tumors are stiffer than the surrounding breast tissues and malignant tumors are much stiffer than benign ones [5–7]. In other words, in vivo identification of the mechanical parameters of normal and abnormal tissues should improve the accuracy of cancer diagnosis. Correspondingly, elastography has been proposed as a method to image the tissues’ elasticity in a quantitative manner. The general basis of elastography is to induce motion within tissue by mechanical stimulation. Conventional medical imaging modalities are then used to measure the spatial deformation, from which the mechanical properties can be extracted. Based on the imaging modalities used, elastography has two major classes: ultrasound elastography (USE) and magnetic resonance elastography (MRE). USE, developed in the 1990s, is the first modulus-imaging modality. It computes the lap between the pre- and postcompression radio frequency ultrasound signals to estimate the tissue’s axial displacement and strain under quasistatic loading [8, 9]. While providing new information for detecting pathological tumors, USE suffers from limited stiffness range as imposed by the minimum resolvable wavelength. The computed image in USE is also restricted by the angular resolution of the transducer and its ability to separate signals from artifacts and noise [9]. Magnetic resonance elastography (MRE) is a second-generation elastography modality that provides higher resolution images and is capable of producing sufficient 3D spatial and contrast resolution [10, 11]. MRE is, however, significantly more costly as a result of the MR imaging procedure and hence is not generally applicable for all patients. From the viewpoint of solid mechanics, the current USE and MRE are insufficient, because both are based on infinitesimal-strain linear elasticity and only very few are capable of considering anisotropic tissue properties. In other words, the large deformation, nonlinear, and anisotropic behaviors of breast tissues (fat and glandular tissues) and tumor have not yet been taken into consideration by USE or MRE. Therefore, the outcomes of USE and MRE may not be sufficiently accurate for the diagnostic purpose.
Motivated by the significance of early detection of breast tumors and the current limitations of mammography and elastography modalities, we have developed a nonlinear elasto-mammography method that takes into consideration of the finite-strain nonlinear properties of breast tissues, in combination with mammography visualization. The development has experienced two stages.
First, a linear elasto-mammography framework was developed to generate the elastograms of breast tissues, by combining the conventional low-dose X-ray mammography with linear elastography framework [12]. Instead of applying ultrasound or magnetic resonance as in the previous elastography research, elasto-mammography uses displacement information extracted from mammography projections before and after breast compression. Incorporating the displacement measurement, an elastography reconstruction algorithm was specifically developed to estimate the elastic moduli of heterogeneous breast tissues. Case studies with numerical breast phantoms showed that the displacement measurement obtained from mammography is sufficient to identify the material parameters of breast tissues and tumors within the framework of linear elasticity.
Then, a nonlinear elastography method was proposed [13]. As discussed above, the current elastography (USE or MRE) reconstruction framework is based on the assumption of linear elasticity theory. The mechanics of biological soft tissues, however, require nonlinear continuum mechanics description [14, 15]. While tissue models based on linear elasticity have been broadly used, they are reliable only when the tissue strain is less than 5% [16], which is much lower than the deformation of soft tissues. Thus, consideration of nonlinearity is essential for elastography in clinical applications. Our development of nonlinear elastography method, for the first time, enables identification of the mechanical properties of soft breast tissues and tumor. To improve the computational efficiency and enhance the stability, a nonlinear adjoint method was introduced. The phantom study demonstrated that the complex nonlinear mechanics of soft breast tissues and tumors can be quantified from 3D displacement and force measured on the surface of the breast.
The objective of the present study is to develop a nonlinear elasto-mammography framework that combines the simplicity of projective X-ray mammography measurement with the accuracy of nonlinear elastography. In Section 2, we present the mathematical derivation, where an adjoint gradient method is modified to consider the projective displacement measurements. Finite-element- (FE-) based numerical simulations are conducted in Section 3 to reconstruct the material parameters of a 3D heterogeneous breast phantom from mammography displacement. Two types of mammography compressive loadings are applied, and the displacements at key points on the tissue interfaces are extracted from mammography projections before and after deformation. In Section 4, the results are presented and the effect of experiential error is investigated.
2. Methods
2.1. Finite-Strain Deformation Equations
Let be a biological object subjected to body force and surface force on boundary . Here, we consider general problems that the body force and surface force are deformation dependent. Following the standard finite-element method, the displacement is discretized as nodal displacement vector , where corresponds to prescribed on and is to be solved from nonlinear equations; that is, on surface (), as described in [13], the FE description of the finite-strain equilibrium equation is The internal nodal force corresponds to the stress of the tissue; that is, it changes with and material parameters but not as it is prescribed. The external nodal force is due to the prescribed surface force and body force in biological object . It changes with the displacement in large deformation. The nodal force is the unknown constraint force on .
A classic quasi-Newton method [17] is employed to solve (1) for . Let be the trial solution of the unknown at the th iterative step. An improved solution can be obtained at the next step, in which is the solution of linear equations: with where the matrices are evaluated at .
2.2. Nonlinear Elasto-Mammography Algorithm
We consider that the biological object is discretized into FE mesh, and the displacement and force are discretized consistently into nodal displacement and nodal force. Experimental measurement for elasto-mammography is displacement. We catalog the measurements as the following. (i) If the force at a node is known, it will be included into which is considered “prescribed” in (1). The corresponding nodal displacement will be considered as unknown in the FE equation (1). (ii) All the other nodal displacements will be in , and the corresponding unknown nodal force will be in . For category (ii), must be considered “prescribed” to fulfill the requirement of the well poseness of a solid mechanics problem.
In our previous elastography method [13], displacements are also measured at some of the nodes associated with and are denoted as . Given material parameters , the unknown displacement and constraint force (which depends on ) will be solved from the FE equation (1). The elastography method thus seeks so that the overall difference between measured and computed is minimum; that is, to minimize objective function: where diagonal weight matrix , with component when the th component of is experimentally measured, or otherwise.
In mammography, however, the measurement of displacement is limited by the projection; that is, only the two components perpendicular to the projection direction are obtainable. Correspondingly, the computed displacement should be projected in the same direction as in mammography and then compared with the mammography measurement . As derived in Appendix A, the projection can be represented by a linear translation of , as , where is a global projection matrix. The objective function for nonlinear elasto-mammography is then
2.3. Nonlinear Adjoint Method
Efficient and robust optimization-based elastography reconstruction schemes request user-supplied gradient . Direct calculation of the gradients involved in the minimization-based parametric identification is difficult, because is an implicated function of . Recently, an adjoint method was introduced to compute the gradient analytically [18–21]. The corresponding nonlinear finite element formulas are shown in Appendix B. Briefly, given a trial will be solved from FE equations (2) and (3), the objective function will calculated by (5), and the material parameters will be updated by large-scale limited memory BFGS (L-BFGS) method with user supplied gradients readily obtained as: where the virtual adjoint displacements and are solved from linear equations: with the tangent stiffness matrix defined in (3). The most significant features of the adjoint method are the analytical formulation, high accuracy, and computational efficiency [22]. Since and its LU factorization have been calculated when solving the FE equation (2), the additional computational expense for in (7) is minimal. Furthermore, it only needs to solve one linear equation (7) regardless of the number of unknown parameters in .
The reconstruction procedure is illustrated in Figure 1. We first establish a numerical FE model of the breast tissue on which external loadings are applied. In order to measure displacement, we compare the mammography projections before and after the deformation. Then, initial guess of the distribution for material parameters is given. Given the external loadings and material parameters, the displacement filed is solved from (1) and is projected to according to the mammography direction. The difference between prediction and measurement are evaluated by the objective function (5). The adjoint field is calculated by (7), and gradients are obtained by (6). The material parameters could be updated by limited-memory BFGS (L-BFGS) optimization subroutine [23]. The iteration continues until a minimization is reached.
3. Numerical Simulations
3.1. Breast Phantom and Forward Problem
We establish a 3D typical breast FE phantom, shown in Figure 2, consisting of the fatty and glandular tissues and a ductal carcinoma (tumor). Boundaries of these regions are described with sets of splines. The mechanical properties of these tissues are described with Fung-type isotropic hyperelastic model [14], whose strain energy function reads where is the Green strain and are material parameters. The parameters are previous determined [13] from ex vivo experimental data of Samani and Plewes [24] as , , ( and are dimensionless, is in kPa) of ductal carcinoma, of fatty tissue, and , , glandular tissue.
Motivated by the breast compression in X-ray mammography, we designed two loadings as detailed in [13]. In the FE model, the base of the breast phantom is fixed. Two paddles are used to apply displacement on the upper surface of the breast. The paddle close to tumor applies tilted compression, and another paddle is fixed to restrict the breast.
3.2. Acquisition Projection Data
For each loading, mammography projections for 3D heterogeneous breast phantom are taken before and after deformation (Figure 2). To mimic the displacement obtainable from mammography, we extract the displacement components in the projection plane at some discrete material points (Figure 3), denoted as . We select three mammography projection directions. With each direction, one projection is made at undeformed state, and one is made at deformed configuration (Figure 2). Then, the displacement components on the projection plane are extracted from a set of landmarks in the tissues by comparison their position in undeformed and deformed projections, as shown in Figures 3 and 4. The landmarks include the top vertex on the upper breast surface (point A in Figure 3), four vertexes of the tumor surface (points B–E in Figure 3), and ten material points on the fat-glandular interface (points A–J in Figure 4). It is noted that the surfaces of tumor and glandular tissue are not smooth so that there are plenty of landmarks that can be used to track the deformation.
To explain the procedure, we use a mammography compression as example. Figure 2 shows mammography projection taken in the same direction with compression applied on the breast. The boundary of the fatty tissue, glandular tissue, and a tumor can be seen in the projection. Then, displacement components on the projection plane can be extracted by comparing the undeformed and deformed projections (Figures 3 and 4). More specifically, the undeformed and deformed projections of fatty tissue and the tumor are registered and shown together for the comparison. The top vertex of fatty tissue, point A, moves to vertex A′ after deformation. Points B–E are vertexes of the tumor in undeformed projection, and they move to vertexes B′–E′ after deformation (Figure 3). On the fat-glandular surface, we select additional ten landmarks that move from A–J to A′–J′, respectively (Figure 4). Thus, by measuring the vector from a point to its deformed position, for example, AA′, the projective displacement components are obtained and recorded as . In addition, it is assumed that there is no slip between the paddles and breast surface during mammography compression. Therefore, the displacement of the material points directly compressed by the paddles is considered known and is added to the measurement .
In summary, we have obtained the following displacement measurements from mammography compression: (i) the top vertex on the upper breast surface and four vertexes of the tumor; (ii) ten nodes on the fat-glandular interface; (iii) material points directly compressed by the paddles. These displacement measurements are denoted as and will be used to identify the material parameters of the tissues.
3.3. Identification of Material Parameters from Displacement Measurements
Having obtained measurement from mammography compression, the inverse problem will be conducted to identify the material parameters of the breast tissues and tumor, with use of an iterative optimization procedure (Figure 1). A homogeneous initial guess of and are dimensionless, is in kPa) is used for all the materials. With a trial , the displacement field is solved from the FE equation (1) and is projected to according to the mammography direction. The difference between prediction and measurement is evaluated by the objective function (5). The gradients are computed with the proposed nonlinear adjoint method. Then, a modified trial will be obtained according to the present and by using L-BFGS minimization subroutine [23]. The iteration continues until a minimization is reached, which corresponds to identified material parameters.
4. Results and Discussion
4.1. Ideal Input
Table 1 shows the initial estimate and reconstructed results, together with the real values for comparison. 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 0.3% ( for tumor) since the effect of the tumor on surface force measurement is the smallest. Reconstructions using different initial estimates have been conducted and very similar results are found, which indicates the efficiency and uniqueness of the proposed nonlinear elasto-mammography using projective measurements. In our study, all numerical experiments reached convergence and had similar convergent profiles. The iteration speed is related with initial estimations. In clinical practice, the initial estimates could be selected based on data of previous patients and experiments. The more reasonable the initial estimates are, the faster the solver got convergence.
In nonlinear elastography [13] and this study, the same nonlinear material model and properties are applied. For ideal input, both frameworks can get convergence and the reconstructed results are very close to the real values. For input with noises, both frameworks could get convergence and have the similar profiles. The parameters in fatty and glandular tissues get convergence faster than these in tumors because the fatty and glandular tissues have bigger impact on surface deformation and measurement.
Convergent loci of the elastic parameters () is plotted in Figure 5. It is observed that elastic parameters of fatty tissue and glandular tissue approach the real values rapidly. After about 50 iteration steps, their relative errors are well within the range of 5%. Then, they experience some minor adjustment. In contrast, elastic parameters of the tumor converge slower. They start to fall to the real values after 300 steps. After 350 steps, all parameters are accurately identified. Reconstructions using different initial estimates have been conducted. Very similar convergent profiles are found, and equally accurate results are obtained. This indicates uniqueness of the proposed elasto-mammography for nonlinear breast tissue properties and efficiency of the reconstruction algorithm.
The slower convergent speed of elastic parameters of the tumor is explained by the roles they play in the deformation due to the applied loadings, as discussed by Liu et al. [18]. In general, parameters with the most significant influence on the deformation are those most easy to identify. The influence of a parameter depends on size and location of the material region it belongs to, as well as characteristics of the deformation. For the present simulations, elastic parameters of fatty tissue and glandular tissue are dominant; those of tumor are much less influential, due to the small size and deep location of the tumor. So parameters of fatty tissue and glandular tissue are more accurately and easily identified than those of the tumor (Figure 5). Therefore, for successful characterization of the tumor, it is critical to apply deformation modes and acquire displacement data that are most affected by the tumor. In this elasto-mammography simulation, displacements of key points on the tumor are extracted from mammography projections, which increase the accuracy and efficiency to reconstruct the elastic parameters, especially for the tumor.
4.2. Multiple Sets of Measurements
Because of the nonuniqueness nature of most inverse problems, it is important to obtain sufficient measurements to reduce the likelihood of nonuniqueness. For 2D isotropic elastography, Barbone and Bamber [25] have shown that one set of displacement and force measurement, especially when measured only on the boundaries, may not provide sufficient information for reconstruction of the distribution of elastic modulus. To enhance the uniqueness of inverse problems, Barbone and Gokhale [26] proposed the feasibility of using multiple displacement fields, and Liu et al. [18] further discussed the use of multiple sets of measurements in 3D anisotropic media. In our previous nonlinear elastography study [13], measurements from four independent titled compression loadings were used to insure stable and unique material parametric reconstruction. In this work, we applied only projective measurements from two breast compression tests and found that the acquired displacement and force data are sufficient for stable parametric reconstruction, even for the small and deeply embedded tumor. This is a significant reduction, as it increases the clinical efficiency, reduces X-ray dose and operation cost, and benefits the patients.
The reduction of necessary loadings is possible because mammography projection provides displacement on the surface of the tumor, which contains direct information of the mechanics of the tumor. Our previous nonlinear elastography study [13] takes only measurement on the breast surface as input. The lack of necessary constitutive information of the tumor in the surface measurement must be compensated by increasing the number of required loadings. In case that the measurement may contain experimental errors, we must use four loadings in the elastography study, instead of two in the present elasto-mammography.
4.3. Iteration Steps
The nonlinear elasto-mammography reconstruction uses an iterative optimization procedure (Figure 1), which is controlled by user-defined criteria. This study employs more strict criteria than in our previous work [13], and it takes about 590 steps to reach the converged reconstruction results. To demonstrate the intermediate results, the uniaxial tensile strain-stress curves of the tumor predicted by the updated material parameters are plotted in Figure 6 at the 1st, 100th, 200th, 300th, and 592nd iterative steps and compared to the real one. It is observed that the reconstructed strain-stress curve approaches the real one rapidly in first 300 iterative steps. After that, the reconstruction only applies some minor adjustment.
In clinical practice, a more tolerable criterion may be applied to control the iterative reconstruction procedure to save computational expense and time. It has been recognized that tissue stiffness plays an important role for diagnosis of breast cancers, as tumors are stiffer than that surrounding breast tissues, and malignant tumors are much stiffer than benign ones [6]. In another word, the stiffness ratio between fatty tissue and tumor, instead of real material parameters, could be used to determine the character of tumors. It is observed in Figure 6 that, starting from the 100th iterative step, the stiffness ratio of tumor to fatty tissue (the lowest curve) increases rapidly, indicating that the predicted mechanical properties of the tumor are well distinguished from the normal tissues for characterizing the tumor. Therefore, from clinical point of view, the iterative reconstruction procedure could be stopped after about 100 steps.
4.4. Input with Noise
The above elasto-mammography reconstructions are conducted using ideal inputs. However, noise is unavoidable in experimental data. To investigate the capability of the proposed nonlinear elasto-mammography modality and algorithm to handle imperfect experimental data, we conduct reconstruction using noisy input, where a randomly selected relative error between ±5% or ±10% is added to each displacement data in . For each noise level, three case studies are conducted. The results are shown as noise 5% (I)–(III) and noise 10% (I)–(III) in Table 1, and the reconstructed tensile strain-stress curves of the tumor are plotted in Figure 7.
It is observed that the strain-stress curves reconstructed with noisy input have similar shape to the ones with ideal input. It is not surprising that curves with 5% noise are closer to the real one than these curves with 10% noise. It demonstrated that, in order to get robust results, we need to make effort to decrease the noise in displacement measurements. It is noted that all the predicted strain-stress curves of tumor, with or without measurement noise, are well distinguished from the curve of fatty tissue (the lowest curve in Figure 7); that is, being much stiffer. That is, even though measurement noise exists, the tumor can be identified by recognizing the difference of stiffness between tumors and the surrounding tissues. This demonstrates that the nonlinear elasto-mammography results are accurate enough for diagnosis of tumors in clinical application.
The previous nonlinear elastography based on surface measurement [13] fails to reconstruct material parameters when ±5% random noise is added to the input. A regularization is required to provide additional constrain. In comparison, the present elasto-mammography yields accurate enough material parameters even with ±10% random noise. The reason is, as mentioned in Sections 4.1 and 4.2, that the displacements extracted on the surface of the tumor from mammography projections contain direct information of the mechanical properties of the tumor, which enhances the robustness of reconstruction and increases the accuracy, in particular of the tumor’s parameters.
4.5. Advantages of Nonlinear Elasto-Mammography
In this study, a nonlinear elasto-mammography framework is developed to incorporate the conventional X-ray mammography for characterization of breast tissue properties. This work extends our previous study linear elasto-mammography [12] and nonlinear elastography [13]. Comparing with previous study, nonlinear elasto-mammography has the following three major advantages.
Imaging techniques: an imaging technique should be selected to measure deformation in elastography. In the proposed nonlinear elasto-mammography, the deformation is measured by conventional X-ray mammography while USE or MRI is applied in nonlinear elastography. Traditional X-ray has advantages of low cost and high resolution, compared with USE and MRI.
Deformation theory: the linear elasto-mammography framework is based on infinitesimal strain deformation theory. However, it is well known that the mechanical behavior of biological soft tissue is nonlinear. In nonlinear elasto-mammography, nonlinear material model and deformation theory are applied so that more accurate results could be obtained.
Inversion techniques: once displacements are measured, an inversion technique is applied to reconstruct elastic properties. In linear elasto-mammography, an adjoint method is applied and then a nonlinear adjoint method is developed for nonlinear elastography. In this study, the nonlinear adjoint method is further improved to enhance the numerical efficiency and stability of reconstruction of elastic properties.
Therefore, the proposed nonlinear elasto-mammography framework has advantage of imaging techniques, deformation theory, and inversion techniques. It combines the simplicity of projective X-ray mammography measurement with the accuracy of nonlinear elastography.
5. Summary
This study presents a nonlinear elasto-mammography method that combines elastography reconstruction and X-ray mammography imaging for the purpose of diagnosis of breast tumors by identification of the finite-strain mechanical parameters of breast tissues and tumors. The displacement information of selected material points is extracted from mammography projections before and after breast compression. Correspondingly, the previously developed nonlinear elastography algorithm has been adjusted with a revised adjoint gradient method to incorporate projection-type displacement measurement. The simulations with heterogeneous breast phantom proved the feasibility of elasto-mammography and tested the efficiency and robustness of the reconstruction algorithm. The simulations show that the deformation of the tumor, depicted by the projected displacement on the surface of the tumor extracted from mammography images, is critical for the success of elasto-mammography reconstruction.
Appendices
A. Displacement Transition between Coordinate Systems
This appendix presents the transition of a displacement vector between global coordinates and projection coordinates. The outcome is the projection matrix in formulas (5) and (7).
To be consistent with computational geometry, we call the projection coordinates as eye coordinates. As illustrated in Figure 8, the global coordinates are denoted as and eye coordinates are . Their direction vectors are and , respectively. As shown in Figure 8, the eye coordinates rotate from global coordinates by three angles: -axis tilt angle , twist angle about eye/original ray , and rotation angle about -axis . It can be shown that
where the rotation matrix is
Now, consider a displacement vector of a material point from undeformed position to deformed position. The global coordinates of are , the eye coordinates are , and their relationship can be derived as
In mammography projection, the displacement component in direction, , is not obtainable, and only and are measured. Therefore, (A3) reduces to
Finally, the FE solution of displacement field , when projected, becomes where is the assemble of according to the FE discretization and assembling methods.
B. Adjoint Method for Gradients of Objective Function
Direct calculation of the gradients of the objective function involved in the minimization-based parametric identification is difficult, because is an implicated function of . An adjoint method will be derived here for efficient and analytical calculation of the gradients. To release the implicit coupling between and , we introduce the constraint (1) into the objective function (5) and obtain a Lagrangian:
where and are arbitrary virtual displacements. In this Lagrangian, and are explicit variables and are no longer coupled. It is noted that and for arbitrary and under the constraint (1). The variation can be expressed as
for which the equality constraint (1) has been applied. Note that the prescribed external force is independent of . Equation (B2) can be further simplified by letting the arbitrary virtual displacement , as
If we select a to let for arbitrary , we obtain a simplest form of , as
Consider that for arbitrary and , we obtain (6) in the text with the following selection of and :
which is (7) in the text.
By introducing the adjoint method, it seems that more equations (B5) and variables ( and ) are involved. But the solution of (B5) is straightforward and the computational cost is minimal, because has been computed and factorized when solving for the displacement as in (3).
The gradients can also be calculated directly as
in which can be computed numerically using finite-different method:
or analytically by solving linear equations:
For finite-strain nonlinear problem, the finite-different method is unaffordable due to the high computational expense to solve (1) for . Solving (B8) is straightforward and is much less expensive for has been computed and factorized. However, (B8) needs to be solved for every material parameters involved; for example, in the exemplar simulations in this work, it needs to be solved nine times because each material has three parameters. In comparison, the proposed adjoint method (B5), (B6) requires only one solution for , regardless of the number of material parameters involved.
Acknowledgment
This work is supported by the US Army’s Breast Cancer Research Program Concept Award (W81XWH-05-1-0461).