A polynomial hyperelastic model for the mixture of fat and glandular tissue in female breast

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN BIOMEDICAL ENGINEERING Int. J. Numer. Meth. Biomed. Engng. 2010; 00:1–?? Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/cnm A polynomial hyperelastic model for the mixture of fat and glandular tissue in female breast Jose. L. Calvo-Gallego*†, Javier Mart´ ınez-Reina and Jaime Dom´ ınguez Department of Mechanical Engineering, School of Superior Engineering, University of Seville, Seville, Spain SUMMARY In the breast of adult women, glandular and fat tissues are intermingled and cannot be clearly distinguished. This work studies if this mixture can be treated as a homogenized tissue. A mechanical model is proposed for the mixture of tissues as a function of the fat content. Different distributions of individual tissues and geometries have been tried to verify the validity of the mixture model. A multiscale modelling approach was applied in a finite element model of a representative volume element (RVE) of tissue, formed by randomly assigning fat or glandular elements to the mesh. Both types of tissues have been assumed as isotropic, quasi-incompresible hyperelastic materials, modelled with a polynomial strain energy function, like the homogenized model. The RVE was subjected to several load cases from which the constants of the polynomial function of the homogenized tissue were fitted in the least squares sense. The results confirm that the fat volume ratio is a key factor in determining the properties of the homogenized tissue, but the spatial distribution of fat is not so important. Finally, a simplified model of a breast was developed to check the validity of the homogenized model in a geometry similar to the actual one. Copyright © 2010 John Wiley & Sons, Ltd. Received . . . KEY WORDS: Breast biomechanics; hyperelastic; finite element analysis; fat tissue; glandular tissue; numerical homogenization 1. INTRODUCTION Breast biomechanics is a field of great interest in the last few years. Procedures as mastectomy, tumorectomy, breast augmentation or reduction have been investigated by many authors, some of whom have developed computational models to study them. For instance, P´ erez del Palomar et al. [1] constructed a finite element (FE) model of the breast from CT images, and applied gravity loads to predict the deformation in supine and standing up positions. Lapuebla-Ferri et al. [2] simulated with a FE model the prosthesis insertion during an augmentation mammoplasty. Roose et al. [3] also studied methods for breast augmentation. They presented a computational model capable of ∗Correspondence to: Jose. L. Calvo-Gallego,Department of Mechanical Engineering, School of Superior Engineering, University of Seville, Seville, Spain †E-mail: joselucalv[email protected] Copyright © 2010 John Wiley & Sons, Ltd. Prepared using cnmauth.cls [Version: 2010/03/27 v2.00] 2J. L. CALVO-GALLEGO ET AL. simulating the postoperative shape of the breast with an accuracy up to 1cm after a subglandular breast implantation. Azar et al. [4] oriented their model to guide clinical breast biopsy. Rajagopal et al. [5] studied how to generate customized FE models by fitting geometrical models to segmented data from MR images. Computational models of the breast must face the problem of the complex distribution of tissues within the organ. The female breast is composed of various kinds of soft tissues, mainly fat, gland, Cooper’s ligaments and skin. It is known that the fibroglandular tissue has a ramified distribution, starting from the nipple [6], but is intermingled with the fat, thus, making the segmentation of both types of tissue from CT or MR images difficult in general [7]. Moreover, the Cooper’s ligaments are difficult to distinguish even with advanced scanning techniques [8]. Therefore, building a computational model of the breast distinguishing all the tissues is a hard task. Another difficulty is the fact that the volume fraction and the spatial distribution of each tissue vary much among women [1]. Because of this, most FE models of the breast do not distinguish between the different tissues present in the breast [1,2,3,4,5]. Instead, they model the breast as composed of only one material (a mixture of all tissues) covered with a thin layer of skin. The question addressed here is if that mixture of tissues is suitable to study the breast, and which are its mechanical properties (as a function of the fat volume ratio) to reproduce the behaviour of the breast with a reasonable accuracy. The material model chosen for the tissues is another key factor, since it has a strong influence on the accuracy of these computational models [9,10]. The breast is composed entirely of soft tissues. These materials usually have a strongly non-linear behaviour because of their internal structure and water content. Their behaviour is usually approximated as quasi-incompressible and hyperelastic. Depending on the importance of dynamic loads, viscoelastic effects can be added [11]. Moreover, depending on the material, they can be considered as isotropic or anisotropic, for example, by modelling the collagen fiber directions [12]. Some information can be found in the literature about modelling the mechanical properties of fat and fibroglandular tissues. All the works model both tissues as isotropic and quasi-incompresible materials. Samani et al. [13] considered them as elastic and determined experimentally the Young’s modulus of samples of both tissues. Some time later, Samani and Plewes [14] determined experimentally their elastic properties, but considering the materials as hyperelastic, with a five terms polynomial strain energy function. Other studies by Samani et al. [15,16] and by O’Hagan and Samani [17,18] have provided the properties of breast tissue with tumoral inclusions. Other authors have tried to obtain the mechanical properties through numerical methods. Usually, they acquired CT or MR images in two body positions, for example, supine and standing up positions. Then, a FE model was constructed and the constants of the strain energy function were determined by fitting the relative displacements between both positions. P´ erez del Palomar et al. [1] fitted the constant of a neo-hookean model for the tissue inside the breast with this technique. This tissue was then covered with a layer of skin, modelled with a polynomial strain energy function and with constants taken from the experimental results obtained by Gambarotta et al.[19]. In the present work Cooper’s ligaments were not modelled. These ligaments can play an important role in the mechanical behaviour of the breast in young women, but as women get older ligaments become slack and its role is less important [20,21]. This work is mainly oriented to simulate the Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm A POLYNOMIAL HYPERELASTIC MODEL FOR THE MIXTURE OF FAT AND GLANDULAR TISSUE IN FEMALE BREAST3 behaviour of the breast in middle-aged to old women. Thus, their effect is not so important, in general. This study is focused on the inner tissue of the breast, which, in a first approximation, is modelled as a homogenized mixture of fat and fibroglandular tissue. The numerical method for homogenization of the mixture is based on the procedure proposed by Temizer and Zohdi [22]. The general objective of this paper is to propose a mechanical model for the behaviour of that homogenized material, providing its properties as a function of the fat volume ratio (termed as f from now on). This is intended to simplify the computational models of the breast and to eliminate the need for segmentation of fat and fibroglandular tissue. To achieve this objective, it is important to study how the distribution of both tissues affects the global behaviour of the breast. Thus, the mechanical behaviour of different mixtures of fat and fibroglandular tissue has been analysed as well as the effect of different proportions and distributions of tissues on the model parameters. 2. METHODS In this section, all the conducted studies and simulations are explained. They have been organised in groups of experiments, with different objectives. For a better understanding a summary table is presented (table I). 2.1. Study A To determine the elastic properties of the homogenized material as a function of the fat volume ratio, a simple FE model has been built using Abaqus FEA®. The model consists in a cube with dimensions 1x1x1 mm3and meshed with 8000 type C3D8H elements (3D eight-noded hexaedral hybrid elements), 20 equally sized elements per edge. The boundary conditions applied to the model depend on the load case considered: uniaxial tension, uniaxial compression or shear, and are explained below. In addition, displacement symmetry boundary conditions have been imposed in all cases in the x= 0,y= 0 and z= 0 faces of the cube (see figure 1). These symmetry conditions would be necessary in a homogeneous material to simulate a uniaxial stress state and were applied to the heterogeneous model for consistency of the homogenization procedure. In tension, a uniform displacement in zdirection has been applied to the z= 1 face of the cube, with no other restrictions. Since the behaviour is assumed incompressible, the deformation gradient tensor enforced is: FT=   1 √λ 1 √λ λ   (1) The stretch, λ, has been increased from 1 to 1.5. In compression, identical boundary conditions were applied, but now decreasing λfrom 1 to 0.7. Finally, in the pure shear load case, uniform displacements were enforced in the yand zdirections, such that the deformation gradient tensor is: Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm 4J. L. CALVO-GALLEGO ET AL. FS=   1 1 λ λ   (2) The stretch, λ, has been increased from 1 to 1.5 in this case. Next, material properties were assigned to each element, corresponding to fat or fibroglandular tissue. This assignment was random, with a uniform spatial distribution, that is, with no bias for any position within the cube, except for accomplishing a given fat proportion (in terms of volume ratio). This random distribution was obtained by numbering the elements of the mesh and assigning the fat properties to a series of elements obtained with the Matlab function rand. This function generates uniformly distributed pseudorandom numbers, which pass the required statistical tests of randomness and independence. The fat and fibroglandular tissue were assumed isotropic, hyperelastic and quasi-incompressible materials, and were modelled with a 5 terms polynomial strain energy function [14]: Ψ = C10(I1−3) + C01(I2−3) + C11(I1−3)(I2−3) + C20(I1−3)2+C02(I2−3)2(3) where I1and I2in equation (3) are the first and second invariant of the right Cauchy-Green tensor, respectively. The constants Cij of the fat and fibroglandular tissues were taken from Samani and Plewes [14] and are shown in table II, corresponding to the fat volume ratios f= 0 and f= 1. Five fat volume ratios were tested: 10%, 30%, 50%, 70% and 90%. Sixteen cases with different random spatial distributions were studied for each fat volume ratio. Figure 1shows one of the 16 cases for each fat volume ratio. Each case was subjected separately to the three load cases (tension, compression and pure shear) by applying the corresponding displacement in Nsmall increments. The evolution of Cauchy stresses during the load (σT zi in tension, σC zi in compression and σS yi and σS zi in pure shear, for i= 1, ..., N) was evaluated with help of the FE solution, as usual in homogenization procedures. The nodal reaction forces in each direction and each face of the cube were retrieved from the FE solution. These nodal reaction forces were summed to obtain the total reaction forces, Rij, in each direction iand each face j, and with them the components of the first Piola-Kirchhoff stress tensor, Pij =Rij/A0, where A0is the area of the corresponding face in the initial configuration. Then, with the well-known relation σ=J−1PFT, the effective Cauchy stress tensor was obtained. It is important to note that the stress states were not uniaxial in the heterogeneous models (or that corresponding to a pure shear state in that particular case). However, the stress components that violated the uniaxiality were negligible (two or three orders of magnitude lower than the uniaxial stress) as it is well-known from homogenization theory. The homogenized tissue was assumed isotropic, hyperelastic and quasi-incompressible, and modelled with the same strain energy function as the fat and fibroglandular tissues (equation 3). The Cauchy stresses of the homogenized tissue can be represented analytically as a function of the constants Cij and the stretch as follows. First, the analytical Cauchy stress tensor is given by: σ=−pI+ 2( ∂Ψ ∂I1 +I1 ∂Ψ ∂I2 )b−2∂Ψ ∂I2 b2(4) Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm A POLYNOMIAL HYPERELASTIC MODEL FOR THE MIXTURE OF FAT AND GLANDULAR TISSUE IN FEMALE BREAST5 where pcan be interpreted as a Lagrange multiplier, identified with a hydrostatic pressure, Iis the unit tensor and b=FFTis the left Cauchy-Green tensor. In uniaxial tension this tensor is: σT=   0 0 0 0 0 0 0 0 σT z an   (5a) The same applies for uniaxial compression, just substituting σC zan for σT zan. In pure shear, this tensor is written in principal axes as: σS=   0 0 0 0σS y an 0 0 0 σS z an   (5b) Calculating the btensors from the deformation gradient tensors of equations (1) and (2), and substituting bin equation (4), the following components of the Cauchy stress tensor are obtained in each case: σT z an = 2C10(λ2 − 1 λ)+2C01(λ− 1 λ2)+6C11(λ3 −λ2 −λ+1 λ+1 λ2 − 1 λ3)+ + 4C20(λ4 −3λ2+λ+3 λ − 2 λ2)+4C02(2λ2 −3λ− 1 λ+3 λ2 − 1 λ4) (6a) σS y an =− 2(−1 + λ2) λ4C11 + 2C20 + (2C02 +C10 −C11 −4C20)λ2+ + (C01 −4C02 −C11 + 2C20)λ4+ (2C02 +C11)λ6 (6b) σS z an =2(−1 + λ2) λ42C02 +C11 + (C01 −4C02 −C11 + 2C20)λ2+ + (2C02 +C10 −C11 −4C20)λ4+ (C11 + 2C20)λ6 (6c) The expression for σC z an is identical to that for σT z an. The subscript “an” stands for analytical solutions. A measure of the error produced in the estimation of the Cauchy stresses with this method can be represented by the squared error, SE, of the three load cases altogether: SE(C10, C01, C11, C20, C02) = N X i=1 (σC zi −σC z an)2+ N X i=1 (σT zi −σT z an)2+ + N X i=1 (σS yi −σS y an)2+ (σS zi −σS z an)2 (7) Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm 6J. L. CALVO-GALLEGO ET AL. where Nis the number of points used in each load case to compare the analytical and the FE solution. The analytical stresses in equation (7) depend on the constants Cij of the model, which can be fitted by minimizing SE in the least squares sense. With this procedure, a set of constants was estimated for each of the 16 random cases. Finally, the average constants of those 16 cases was obtained for each fat volume ratio. 2.2. Study B After estimating the constants of the homogenized material for each fat volume ratio, other material distributions and geometries have been simulated to test the behaviour of the homogenized material. In study B, the same procedure followed in study A has been conducted, but meshing the cube with less elements. By reducing the number of elements, it is more likely to have larger pieces of a single tissue. Such distribution would be closer to reality if both types of tissue were not as intermingled as in the model of study A. Then, the aim of the study B was to analyse the influence of these larger pieces of a single tissue in the overall behaviour of the homogenized material. Only the case of 50% fat volume ratio has been analysed. Two models, with 1000 and 125 elements, respectively, have been tested and compared to the original model of 8000 elements. The selection of these element sizes (between 50 and 200 µm) is justified since a terminal duct lobular unit (TDLU) has an approximate diameter of 200-600 µm and the little lactiferous ducts, which join those TDLUs, have an approximate diameter of 40-60 µm [23,24]. During breastfeeding, these ducts increase their size, but that can be considered a marginal state of the female breast. 2.3. Study C In study C, a simplified geometry of the tissue inside the breast, consisting in a spherical cap of 8cm in radius and a height of 6cm was built and meshed using 93312 type C3D8H elements. Two models were analysed for this geometry, one distinguishing between fat and fibroglandular tissue, named ramified model, and the other one without differentiation, i.e. using the homogenized material properties. The ramified model is shown in figure 2and has 45.3% of fat. To construct this ramification, the point N1 (see figure 2) was established as the origin of a number of beams equally spaced over the spherical cap. The glandular tissue was assigned to the elements within the beams, whose width was arbitrarily set, such that gland and fat were approximately divided in equal proportions†. This distribution tries to mimic, in a simplified manner, the ramified distribution of glandular tissue in an actual breast, in which the gland is more abundant near the nipple, where the lactiferous ducts converge. It must be clear that the intention of this study was not modelling an actual breast, in which case the skin (not included in this model) should be added and more complex boundaries conditions should be considered to model the interaction with the surrounding tissues and organs. The aim of study C was to compare the behaviour of both models (homogenized and ramified) under gravity loads, so to check the validity of the homogenization procedure when applied in a geometry similar to the actual one. †The 50-50% is the most compromising proportion from the homogenization point of view. This homogenization procedure is more accurate as one constituent predominates and is exact for a 0-100% proportion, obviously. Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm A POLYNOMIAL HYPERELASTIC MODEL FOR THE MIXTURE OF FAT AND GLANDULAR TISSUE IN FEMALE BREAST7 As boundary conditions, all the displacements were constrained in the back plane of the spherical cap, i.e. in the contact with the chest wall, like in other studies [1,2]. Both models were subjected to gravity loads simulating three body positions: supine (gravity acting in the negative xdirection of figure 2), prone (gravity acting in the positive xdirection) and standing up (gravity acting in the negative zdirection). The displacements of the nipple (point N1) in both models (ramified and homogenized) were compared for each load case. The densities considered for fibroglandular and fat tissues were 1.04 g/cm3and 0.93 g/cm3, respectively [25,26,27]. For the homogenized material, the density was calculated using the rule of mixtures. 2.4. Study D Finally, in study D, the same procedure followed in study C (with the same models and materials) was conducted, but now changing the boundary conditions. In study C, all the displacements were constrained in the back plane of the spherical cap. In study D, the displacements were constrained in half of the back plane (in figure 2, if a line is drawn from N2 in z direction, the nodes of the back plane which are to the left), while in the other half the nodes were connected with a certain stiffness to a rigid wall (the same stiffness in directions x, y and z). It is known that the connection between the breast and the chest is carried out through muscles and ligaments, and they have some flexibility. This precise configuration of dividing the breast in two halves was selected for it produced a deformed shape that resembled the actual one more closely than that of the study C. To some extent, the rigid connection of the left part is mimicking the constraining effect of the sternum. The results shown here correspond to a stiffness of those connectors of 0.05N/m, which was adjusted to produce a displacement of the nipple similar to the actual one. Anyhow, this study did not pursue a precise modelling of the boundary conditions of the breast, but studying the behaviour of the bulk material of the homogenized model in a more realistic scenario. 3. RESULTS 3.1. Study A The constants of the polynomial strain energy function were fitted for each of the 16 cases studied per fat volume ratio, by minimizing SE given in equation (7). The mean and standard deviation for these constants were calculated for each fat volume ratio and are shown in table II. The root-mean- squared errors (RMSE) for the 16 cases are also shown in table II. Linear or quadratic regressions were selectively proposed between the 5 constants and f(see figures 3and 4). 3.2. Study B A comparison of the random distribution models for different number of elements: 8000, 1000 and 125; was made for f= 0.5. The constants of the polynomial model obtained are shown in table III. Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm 8J. L. CALVO-GALLEGO ET AL. 3.3. Study C The results of the spherical cap geometry are presented here. The displacement of the nipple (node N1, see figure 2) in both models (ramified and homogenized) for the three body position simulated are compared in table IV. The relative difference between those displacements is also shown in the table. 3.4. Study D The results of the spherical cap geometry, with more flexible boundary conditions, are presented here. The displacement of the nipple in both models for the three body position simulated are compared in table IV. The relative difference between those displacements is also shown in the table. 4. DISCUSSION The small errors found in the least squares fitting of the constants with the three load cases applied (see table II) show that the polynomial model is suitable to reproduce the behaviour of the homogenized tissue, provided that the individual compounds are well modelled with a polynomial model with the constants given by Samani and Plewes [14]. The mean values of the constants vary much with the fat content, confirming that the fat volume ratio is a key factor in determining the properties of the composite tissue. In addition, the standard deviations are small, indicating that the distribution of fat is not so important, at least in the random distribution. A linear regression between the constants and the fat volume ratio fits quite well, except for C01 and C11, for which a quadratic regression fits better. The influence of the mesh size in the homogenization procedure was analysed in study B. The constants of the strain energy function were very similar in all cases, confirming that the homogenization is also valid in a breast in which larger pieces of a single tissue exists. In homogenization techniques, the application of displacement boundary conditions to obtain an average stress provides an upper bound of the effective stiffness. Alternatively, the application of traction boundary conditions to obtain an average displacement provides the lower bound. For this reason, three alternative loads were applied: constant tension, compression and pure shear stress in the corresponding faces of one model with 50% of fat proportion. The stretches, λ, were computed in this case, resulting in negligible differences of less than 0.03%. In the spherical cap (study C) the relative differences in displacement between both models (ramified and homogeneous), shown in table IV, are very small, thus confirming the validity of the homogenized model. Despite these good results, the magnitude of the displacements seems low compared to actual displacements in a real breast [28]. This is likely due to the simplistic boundary conditions imposed in the model (fixing the nodes in contact with the chest wall). The actual connection of the tissues is likely more flexible and makes those displacements to be larger. That is the reason why a second model with different boundary conditions (study D) was analysed. As mentioned before, the intention was not modelling very precisely the actual boundary conditions of the breast, but checking if the behaviour of the bulk homogenized tissue is still valid with boundary Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm A POLYNOMIAL HYPERELASTIC MODEL FOR THE MIXTURE OF FAT AND GLANDULAR TISSUE IN FEMALE BREAST9 conditions that produce a deformation state closer to the actual one. And so it is, given that in study D the displacements obtained with the homogenized model are very similar to those obtained with the ramified model. The usefulness of this paper is linked to the availability of assessing the fat volume ratio within the breast in a straightforward way, at least approximately. Several techniques exist to measure that fat volume ratio. Visual assessment through mammography images is probably the simplest one, but it can only provide a qualitative measurement. With this technique an individual can be sorted out into categories, e.g. 0-25%, 25-50%, 50-75%, 75-100%, [29]. The segmentation’s accuracy can be easily increased if some mathematical algorithms are applied [30,29]. With the use of MR imaging and tissue segmentation the accuracy can be increased further, though paying a high price in manipulation and computation time [31]. For some purposes, the four aforementioned categories can be enough in view of the variation of the material constants with the fat content. Therefore, a quick visual assessment of the fat volume ratio through mammography and the use of the proposed homogenized model can strongly simplify the problem of patient specific breast models. Finally, the limitations of the study are presented. The interfaces between fat and gland have been modelled as a rigid union, i.e. the contact between both materials has not been taken into account. Moreover, the model depends on the suitability of the properties of fat and gland given by Samani and Plewes [14], which have been used to develop the present model. It is important to note also that the strain energy function used by these authors is not polyconvex. However, this fact does not necessarily mean that this function should be ruled out, as stated by Hartmann [32,33]. In particular, this author showed that the only requirement for polynomial models to produce monotonous stressstrain curves (that is, with a physical meaning) is that the constants are positive, like they are in the homogenized model for any fat content. Finally, the homogenization procedure carried out here provides a material model which is adequate to simulate the bulk behaviour of the mixture of gland and fat, but it would fail to describe the behaviour under very localized loads, like a needle insertion. In that case, the size of the needle could be of the same order of magnitude as that of the breast constituents. The structural organization of the tissues might play an important role and a microscopic model would be needed in such case. 5. CONCLUSIONS A mechanical model has been proposed for the mixture of the female breast tissues as a function of the fat volume ratio. A numerical homogenization technique has been used to fit the constants of a polynomial hyperelastic model. Different random distributions of fibroglandular and fat tissues have been tested providing consistent results. The elastic response is highly dependent on the fat volume ratio, but it does not depend on the particular distribution of fat and gland. This conclusion has been confirmed in a model with a similar geometry to a real breast, where the fibroglandular tissue has been arranged in a ramified manner. The proposed homogenized model is suitable and constitutes a strong simplification of the problem. Certainly, it reduces the cost of modelling the breast, since no segmentation of fat and fibroglandular tissue would be necessary anymore. Instead, apart from the skin, not considered in this study, only the soft tissues inside the skin must be identified and modelled as a whole and the fat volume ratio measured or estimated. This simplification is important Copyright © 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm