On the three-dimensional mechanical behavior of human breast tissue

2.1. Overall approach

Our goal herein is the elucidation of three dimensional (3D) mechanical properties of human breast tissue and development of a novel constitutive model for its unique behaviors. We utilized our full 3D numerical-experimental approach [Avazmohammadi et al., 2018] to provide the necessary experimental data and parameter estimation pipeline. We accomplished this by first acquiring human cuboidal breast tissue samples from a cohort of breast cancer patients, making sure to take samples from healthy portions of the breast away from tumor sites. These samples were then subjected to several deformation protocols using a triaxial testing system. Based on initial observations of the resulting data, we formulated a constitutive model, and utilized finite element analysis optimization to determine the optimal values of the model parameters. We also verified the predictive capability, based on a subset of the measured deformation protocols. This fitted model was then analysed to determine its degree of anisotropy which we can then use to infer the inherent anisotropy of the tissue itself.

2.2. Tissue sourcing and anatomical selection

Normal human breast tissue samples were collected from a cohort of breast cancer patients undergoing breast surgery, in compliance with the institutional review board at The University of Texas MD Anderson Cancer Center (PA16-0364). Anatomically oriented 2×2×2-cm samples were obtained from multiple healthy regions in the breast in order to bracket spatial variations in mechanical properties and classified into “Adipose” and “Fibroglandular” groups by our pathologist depending on their relation to physical breast anatomy (n=3 for each group). Specimens were first inert dye marked to identify the anatomic axes (Figure 1-A). Before testing, each specimen was further trimmed to 1×1×1-cm cuboids using microtome blades, with the edges of the specimens aligned to the medial-lateral (Med-Lat), superior-inferior (Sup-Inf), and anterior-posterior (Ant-Pos) directions of the breast (Figure 1-B).

2.3. Experimental setup

To measure the full 3D mechanical response of the tissue, we used a previously developed triaxial testing technique [Avazmohammadi et al., 2018]. Briefly, specimens were mounted by attaching nine-pin arrays to each face of the cuboid with cyanoacrylate gel-based adhesive (Figure 1-C)). Load cells within the attachment systems of the device were used to record the 3D force response of the tissue in real time as deformations were applied, controlled with a custom LabVIEW virtual instrument (National Instruments, TX, USA). Due to the pin attachments of our samples, the deformation of each sample deviated from the homogeneous deformation we sought to impose, since, while the displacement boundary conditions on the pins are consistent with these homogeneous deformations, the zero traction boundary conditions away from the pins allow each specimen side to deform relatively freely. Therefore, these deformations must be simulated numerically rather then assessed analytically. This approach allows for a simple specimen attachment setup but more importantly complete accounting for any fibrous specimen heterogeneity Li et al. [2020b,c]. For simplicity, we report the homogeneous deformation gradient consistent with the imposed pin displacements, with these being identified with the X or 1, 2 or Y, and 3 or Z directions respectively for orientation of the specimen within the mechanical testing device (Figure 1-C)). All testing was performed with the specimens immersed in phosphate buffered saline maintained at 37°C.

2.4. Deformation protocols

Bearing in mind the use of an inverse modeling approach required and the relatively unknown properties of the human breast tissue, deformation protocols were chosen to simultaneously explore the multi-axial behavior of the tissue and produce the requisite data for robust parameter estimation (Table 1). Specifically, the tissue’s mechanical response was assessed in both two-axis tension-compression deformations in which one axis is extended while another axis is simultaneously compressed such that the unconfined third axis remains unstretched and single-axis (unconfined) compression (Figure 2). This is known as a ‘pure shear’ mode. Additionally, these loading paths were imposed on varying combinations of specimen axes to differentiate the tissue’s response with respect to anatomical orientation. We further note that use of simultaneous states of tension and compression allowed determination of bimodular material behavior. In addition, we conducted simple compression tests along each anatomical axis (Figure 2). This comprehensive multiaxial evaluation was necessary, since previous studies focused primarily on the tensile behavior of human breast tissue, despite the central role of compression in mammography and elastography techniques. In practice, each protocol was repeated over 10 cycles to precondition the tissue.

2.5. Data post-processing

We note that while some variability in the experimental data occurred, we observed a high degree of consistency in the same anatomical location sub-group, with a standard error of ≤ 20% at peak deformation. As our main focus was on determining an appropriate form and reasonable estimate of the material parameters, we averaged the data from each individual protocol for each sub-group to create tissue group averaged response. This approach had the added benefit of reducing the effect of measurement noise a priori, and has been successfully used by our group in previous soft tissue studies to develop accurate group responses Sacks and Chuong [1998].

We note too that due to Central Limit Theorem and assuming that the model parameters are independent random variables, uncertainty in the averaged data can be well represented in form of Gaussian distributions, even if the uncertainty in the original data cannot be. The functional independence of the different parameters appearing in our model is established in Appendix C. Solving the inverse problem using the average of the experimental results for each test protocol thus leads to a maximum a posterior estimation of the final distribution. The study of the input uncertainties on the final distribution is very important but is out of scope of the current work, as obtaining concrete estimates of the distribution of physical properties across the population would require a greater number of specimens. This averaged data $\widehat{\mathbf{S}}$ was then used as the basis for all parameter fitting and validation.

2.6. Constitutive model

Because we sought to model any observed inherent anisotropy that may be present in human breast tissues, we decided to choose an appropriately general model form that would allow simulate a range of fit isotropic to fully anisotropic behaviors, while also facilitating quantification of the degree of observed anisotropy. As is commonly done in the modeling of soft tissues, we assumed that the size of the specimen ensured the mechanical and structural properties were approximately homogeneous at the length scale of study. The deformations were represented by a map $\mathbf{x}=\phi (\mathbf{X})$, where $\mathbf{x}$ parameterizes the loaded configuration, and $\mathbf{X}$ parameterizes the unloaded configuration. The breast was modeled as a (pseudo)-hyperelastic material Fung [2013] whose deformation can be locally characterized by the deformation gradient tensor $\mathbf{F}=\text{Grad}\mathbf{x}$. We represented the strain energy in the tissue $\psi (\mathbf{E})$ as a function of the Green-Lagrange strain tensor $\mathbf{E}=\frac{1}{2}({\mathbf{F}}^{\mathrm{T}}\mathbf{F}-\mathbf{I})$, which induces a second Piola stress tensor $\mathbf{S}=\frac{\partial \psi }{\partial \mathbf{E}}$. is common with biological materials, we assume incompressibility because of the high water content present in the tissue, requiring a Lagrange multiplier that augments the stress

where $J-1=\text{Det}\mathbf{F}-1=\sqrt{\text{Det}(2\mathbf{E}+\mathbf{I})}-1=0$ is the incompressibility constraint. For computational purposes, we replace the Lagrange multiplier term with an energy penalty term of the form $K(J-1{)}^2$. Specifically, we adopt a modified anisotropic nonlinear neo-Hookean strain energy density of the form

where $\text{Tr}(\mathbf{E})=\frac{1}{2}(I_1-3)$, and $\alpha =\{\mu ,ℳ,K,a,b\}$ are the material model parameters, $\mu$ being a neo-Hookean like shear modulus in the small strain limit, $ℳ$ being a fourth order tensor of material parameters, $K$ being a bulk modulus, and $a$ and $b$ determining the stronger nonlinearity with respect to straining in the Ant-Pos ($x_3$) direction. This form was inspired by early Fung models, which were the sum of a quadratic form and the exponential of another quadratic form [Fung, 2013]. We add a linear neohookean term and an incompressibility constraint to model possible tension-compression asymmetry, and restrict the two quadratic forms to reflect our observed anisotropy. To simulate near incompressibility, we let $\frac{K}{\mu }\gg 1$ so that the volumetric penalty is the dominant contribution to the energy, even though the bulk modulus is still a parameter that we are fitting. This form was chosen based on a preliminary examination of the experimental fibroglandular data, which revealed a strong nonlinearity in the Ant-Pos direction (here identified with the $x_3$ coordinate) that was not present in the other directions, together with otherwise mild anisotropy. For the adipose tissue, the exponential term did not improve fits, and accordingly, we take $a=0$ for adipose tissue, which also eliminates $b$ from the model.

Additionally, in choosing the form of $ℳ$ we utilized the hypothesis that the breast tissue is orthotropic with respect to the anatomical axes, as we expect that any anisotropy present in human breast tissue most strongly present with respect to these axes. Therefore, by restricting the model’s form to orthotropy with respect to these axes, we captured the dominant portion of anatomical anisotropy. A secondary advantage of this is that we reduce the number of material parameters in our model, which both prevents overfitting and reduces the computational cost of optimizing parameters.

Explicitly in Voigt notation, the Green’s strain tensor can be written as

after which, the tensor of material parameters can be written as

This term, being the quadratic coefficient of the strain energy, can be directly identified with the usual tensor of elastic moduli in linear elasticity, and hence can be decomposed into an isotropic component, and a purely anisotropic component. This allows us to quantify the degree of low order anisotropy by considering the fraction of the norm of the tensor of elastic moduli contributed by the anisotropic component. We note additionally that when making this identification, the term $M_{33}$ must be augmented with the addition of $ab$ and a volumetric term involving $K$ to account for the low order contribution of the exponential term and incompressibility penalty, since the expansion of the energy to quadratic order is equal to

Finally we note that the presence of a linear term in the strain energy density would suggest the presence of residual stress in the reference configuration, since this term directly corresponds to stress in the unloaded configuration. The isotropy of the linear term results in a residual stress that is a pure pressure. This means that the residual stress that would be predicted is counteracted by the incompressibility constraint, and hence this linear term does not create residual stress in the reference configuration. This term does however, create a tension-compression asymmetry as demonstrated in Appendix B, since it is an odd function of the Lagrange strain, eliminating the need to employ a bimodular theory as was done in Klisch [2007].

2.7. The Forward Problem

Due to the nature of the inverse problem that necessitates many forward solves plus the complex nonlinearity in the model finding an efficient solver is crucial. Toward this aim, a Newton-Krylov solver preconditioned with Algebraic Multigrid (AMG) was used with GMRES kernel update every 30 iterations. For smaller deformations a relaxation factor equal to 1 worked well, but for larger deformations reducing the relaxation factor to 0.5 was necessary to handle the increasing nonlinearities. Our experience showed non-preconditioned Krylov solvers were incapable of solving the linear problem, and direct solvers were too expensive to compute as expected. Moreover, an advantage of using iterative solvers for optimization problems is the use of the strong Wolfe condition which allows adaptive updating of the convergence criteria while the optimizer gets closer to the target. In light of the boundary conditions resulting from the nine-pin attachments in mechanical testing, we employed a finite element (FE) representation of the tissue specimen mounted to the triaxial testing device, featuring a volume mesh of tetrahedral elements with surface node subsets on each cuboid face to replicate the nine-pin attachments. Details of this approach are given in Avazmohammadi et al. [2018].

2.8. Inverse modeling, parameter estimation, and validation

The optimal material parameters $\alpha$ for both constitutive model forms (i.e. including or excluding the exponential term) were computed with an inverse modeling framework based on our previous approaches [Avazmohammadi et al., 2018, Li et al., 2020a]. A subset of the tri-axial deformations containing the pure shear protocols and one simple compression protocol was simulated in parallel forward simulations, with each loading path in the set prescribed as Dirichlet boundary conditions to the attachments. The modeled stress $\mathbf{S}(\alpha )$ was computed for a given set of material parameters and compared to the experimentally measured data $\widehat{\mathbf{S}}$ with the cost function (Eq. 6)

where we sum over the indices $m,n$ and $q$ corresponding to the deformation modes, relevant cuboid faces, and experimental data points, and $\Vert \cdot \Vert$ is the ${\mathrm{L}}^2$ norm. The optimal parameters ${\alpha }^{\ast }$ were estimated by minimizing $\Phi$ (Eq. 7):

with the reflective trust region method. Gradients were computed by the central finite difference scheme fed with complex numbers. By using complex numbers both the initial value of the governing PDE and its perturbed value can be computed at a point with a single evaluation which nearly (though not exactly, due to the additional overhead associated with using complex instead of floating point numbers) halves the computational cost. The trust-region size was updated iteratively as needed, especially when the optimizer was far from the optimal parameters. A simplified single element analytical solution was used calculate an initial guess for the optimal parameters, a critical step in the convergence of any Newton-based method. As an extension to our previous approaches, we converted the entire finite element simulations to the open source computing platform FEniCS (fenicsproject.org). We used the Stampede system at the Texas Advanced Computing Center. Once the set of optimal parameters ${\alpha }^{\ast }$ is determined, the remaining triaxial deformations are simulated, and the predicted results are compared with the experimental data, to demonstrate the predictive capacity of our model.

3.1. General observations

The stretch stress-stretch responses from all protocols applied to adipose and fibroglandular breast tissue specimens exhibited mild nonlinearity across all deformation modes. Beyond 20-30% deformation, the initially soft and mostly linear curves of protocols 1 – 6 showed a highly nonlinear increase, consistent with other work on adipose tissue mechanical response. Examining the data from multiple cycles showed that preconditioning effects were most prevalent in the first cycle, whereas subsequent cycles had a more reproducible response with less hysteresis.

Although the response of the adipose tissue was only mildly anisotropic, mechanical testing results revealed stronger anisotropy in the fibroglandular specimens, as can be observed from the disparate high-strain stresses obtained protocol 6 in comparison to protocols 4 – 5. In particular, deformations along the Ant-Pos direction showed greater maximum stress compared to other directions. The fitted parameters for both models are presented in Table 2. Additionally, these parameters were able to fit both compressive and tension behavior for both tissues, indicating that a bimodular model was unnecessary.

3.2. Adipose tissue

The mildly anisotropic model was able to achieve a reasonable fit to the mechanical data of adipose tissue specimens (Fig 3a). The fitted model predicts the stress results of the deformation modes not included in the fitting data reasonably well (Fig 3b).

3.3. Fibroglandular tissue

The overall mechanical response of fibroglandular tissue differed from that of the adipose tissue specimens, achieving greater maximum stresses in tension-compression modes but generally lower maximum stresses in compression. There was also a more pronounced anisotropic mechanical profile across protocols, compared to adipose tissue. The addition of the anisotropic exponential term allowed our model to fit the data well (Fig 4a). Additionally, this fitted model was able to predict the stress in deformation modes outside of the fitting data (Fig 4b).

3.4. Compressive behavior of breast tissue

The form of the fitted model indicates that fibroglandular tissue exhibits a strong stress-strain nonlinearity with respect to simple compression in the Ant-Pos direction in comparison to either the Med-Lat or Sup-Inf directions (Fig 5). This initially appears at odds with the common convention regarding fibrous tissues, which is that fibers are incapable of bearing compressive loads. However, the incompressibility constraint allows us to recast this term as a function of the positive strains in the Med-Lat and Sup-Inf direction. Hence if this nonlinearity is due to the fibrous content, the form of the data suggests that these fibers lie in the plane spanned by the Med-Lat and Sup-Inf directions.

3.5. Anisotropy

We seek to quantify if there are differences in the anisotropy present in adipose compared to fibroglandular tissues. We do this by examining the stiffness matrix present in both models, and determine the relative degree of anisotropy present in these matrices. For comparison, the most general isotropic quadratic function of the Green Lagrange strain has a coefficient matrix that takes the form

Stiffness matrices themselves span a vector space, with the matrices corresponding to isotropic materials spanning a two-dimensional subspace of this vector space. We want to decompose the total space into a purely isotropic subspace, and its orthogonal complement, which we refer to as the anisotropic subspace. We can then project our optimized stiffness matrices onto these two subspaces and compare the relative magnitude of their anisotropic components. The definition of the orthogonal complement and the projection onto these subspaces was accomplished using the natural inner product induced by the referential metric tensor, which differs from the Frobenius norm on 6x6 matrices that one might naively employ.

With this norm, we take the quadratic component of our fitted strain energy, and compute its projection onto the two-dimensional isotropic subspace, and the anisotropic subspace in a way that agrees with the intrinsic geometry of our problem. We then express our model’s quadratic coefficient tensor as the sum of an isotropic term and an anisotropic term. This allowed us determine the degree of our model’s anisotropy at quadratic order by computing the magnitude of the anisotropic component and comparing it to the total norm of the coefficient tensor. Performing this analysis, we determined that our fitted fibroglandular model’s relative anisotropy is 2.04 times the relative anisotropy of the adipose model, providing a quantitative measure of how much more anisotropic fibroglandular tissue is than adipose tissue.

Additionally, we note that the vast majority of each model’s coefficient tensor (over 99%) lies in the subspace parameterized by the bulk modulus, indicating that the assumption of incompressibility for both of these tissues is well justified. Removing this component, we find that the overwhelming majority of both models’ remaining strain energy is anisotropic (99.9% for adipose tissue, and 95.8% for fibroglandular). One might expect that these deviatoric anisotropies would mimic the pattern we see in the total anisotropies, in that we would expect that the fibroglandular tissue’s deviatoric anisotropy would be higher than that of the adipose tissue, however here we observe a reversal. This indicates that while the isotropic component of both tissues is primarily due to incompressibility, the fibroglandular tissue has a comparatively larger isotropic stiffness in shear, relative to the adipose tissue. Therefore, when the volumetric component is totally removed, the fibroglandular tissue’s energy possesses a comparatively larger remaining isotropic component, and hence a smaller deviatoric anisotropy. Furthermore, we found that the fibroglandular tissue displayed stronger anisotropy at higher strains than the adipose tissue, which we modeled by the imposition of the exponential term. While did not rigorously quantified the anisotropy of this term as we have with the rest of the parameters in the model since it does not lie in a clearly defined finite dimensional vector space, we note that it is transversely isotropic with the Ant-Pos direction as its distinguished direction. Overall, these results indicate that beyond the incompressibility constraint, which is isotropic, these tissues display complicated anisotropic behavior that cannot be captured by an isotropic model.

4.1. General findings

In the course of this study, we sought to examine and quantify the mechanical behavior of human breast tissue, and in particular, we sought to assess and quantify the anisotropy present in this behavior. Additionally, we specifically wanted to capture breast tissue’s behavior in compression along the anatomical axes, since most previous studies only studied breast tissue’s behavior in tension. We did this by taking samples classified as either “adipose” or “fibroglandular,” imposed multiaxial deformations upon these, and then fit an anisotropic model to each of these collected data sets. Based on these fitted models, we quantified the relative incompressibility and anisotropy of these tissue types, and found that these tissues can be modeled as incompressible, and that their anisotropy could not be neglected. Not only did we find that both tissue types were anisotropic, we found that fibroglandular tissue displayed substantially stronger anisotropy, necessitating an additional exponential term. Despite this feature, a bimodular model was unnecessary as both compressive and tensile components of stress were well predicted by both models. These results suggest that properly identifying the orientation of strong anisotropy in breast tissue arising from its microstructure significantly improves robustness of its characterization.

We also found qualitative differences in the mechanical behavior of our fitted models. Under the imposition of incompressibility, our fitted adipose model’s energy has a minimum at zero strain, while our fitted fibroglandular model’s energy has a saddle point at zero strain, with local minima nearby. Outside of a small region around zero strain, both of our models’ fitted energies are convex. This allows our fibroglandular model to behave incredibly softly at low strains, matching our experimental data, while providing stability at larger strains. While this lack of convexity at zero strain may initially seem concerning, the difference between the energy at zero strain and the minimum energy is very slight, which yields a model that is capable of capturing the softness of fibroglandular tissue at low strain. Furthermore, because the region of fibroglandular nonconvexity is small, we retain the advantages of a convex model at finite strains; our adipose model is convex so we don’t encounter a similar issue with adipose tissue. This analysis is presented in Appendix A.

4.2. Links to breast tissue structure and organ level function

Our results suggest that human breast tissue is complex in its mechanical response, exhibiting both anisotropy to various degrees and heterogeneity. Additionally, linking this anisotropy to the physical structure of the tissue is necessary; validation with histological analysis of the specimens remains to be performed. Furthermore, the mechanical and structural differences between adipose and fibroglandular regions have yet to be more thoroughly assessed, and more rigorous comparison of the constitutive models have yet to be performed with metrics like the Akaike information criterion. This suggests that further studies more thoroughly mapping this heterogeneity and anisotropy, as well as quantifying the variance in material properties and structure across individuals are necessary. Additionally, a more detailed characterization of the in vivo mechanical environment of breast tissue can serve to guide these studies more efficiently. By performing these studies and having a thorough understanding of the complex mechanical behavior of human breast tissue, both at the general population level and the patient-specific level, we can accurately predict healthy tissue’s behavior during deformation-based diagnostic techniques, and hence identify pathological behavior. This information will also improve confidence in the accuracy of simulation data, such as internal stresses and strains, when direct measurement is impossible or invasive. Knowledge of these internal stresses and strains experienced by the breast can be used to guide the development and systematic improvement of tissue reconstruction techniques, synthetic implants, and surgical interventions.

4.3. Summary

In this study, we developed the first specialized 3D model of the mechanical behavior of human breast tissue using comprehensive 3D numerical-experimental approach. This novel model of adipose and fibroglandular specimens revealed region-dependent differences in tissue mechanical properties, which supports the need for incorporating heterogeneity in whole-breast simulations. Additionally, the form of this model allowed us to quantify the degree of anisotropy in both tissues, allowing us to precisely compare the qualitative differences in the fitted models. Despite this, we did not find it necessary to employ a bimodular model to capture the behavior of either adipose or fibroglandular tissue. Our fitted parameters for both models suggests that both adipose and fibroglandular tissues are effectively incompressible, and so incompressible models are well justified. The fitting of anisotropic models also indicates that increasing fibrous content in the breast tissue could indeed drive mechanical anisotropy at higher strains; however, for low fibrous content and at low tissue strains, the anisotropy of the tissue is well described as a low order phenomenon. We found that this anisotropy in all cases constitutes the majority of human breast tissue’s isochoric response, hence perhaps the primary finding of our work is the dominance of breast anisotropy in volume-preserving motion. The findings presented in this work serve as a critical first step toward physiologically accurate 3D modeling of the breast and highlight important considerations for the constitutive modeling of soft tissues.