Experimental characterisation of breast tissues and its application to a numerical model of a healthy breast

EXPERIMENTAL CHARACTERISATION OF BREAST TISSUES AND ITS APPLICATION TO A NUMERICAL MODEL OF A HEALTHY BREAST by JOSE LUIS CALVO GALLEGO A thesis submitted in fulfilment of the requirements for the degree of DOCTOR OF PHILOSOPHY presented at the ENGINEERING SCHOOL OF THE UNIVERSITY OF SEVILLE supervised by the professors JAVIER MART´ INEZ REINA and JAIME DOM´ INGUEZ December 2016 Nomenclature Throughout this text, scalars, vectors and tensors are going to be used. For scalars, lowercase Greek and uppercase and lowercase italics letters are going to be used (i.e. t,α,J). For vectors referring to the material configuration, uppercase bold-face Latin letters (i.e. N,X) and for vectors referring to the current configuration, lowercase bold-face Latin letters (i.e. n,x). For second-order tensors in the reference configuration, uppercase bold-face Latin letters (i.e. S,E) and for second-order tensors in the spatial configuration, lowercase bold-face Greek or Latin letters (i.e. σ,e). For third-order tensors, uppercase boldface calligraphic letters (i.e. E). For fourth-order tensors in the material configuration, uppercase blackboard Latin letters (i.e. C ) and for fourth order tensors in the spatial configuration, lowercase blackboard Latin letters (i.e. c ). Symbol Description β∞ jStrain energy factor ΓjTensor which determines the relaxation/creep behaviour λaPrincipal stretches, with a= 1,2,3 ¯ λaModified principal stretches, with a= 1,2,3 ρ0Density in the reference configuration ρDensity in the current configuration ¯σRecorded experimental Cauchy stress ˜σTreated experimental Cauchy stress σCauchy stress tensor σvol Volumetric part of the Cauchy stress tensor σiso Isochoric part of the Cauchy stress tensor ¯σFictitious Cauchy stress tensor τiRelaxation time constants τKirchhoff stress tensor τiso Isochoric part of the Kirchhoff stress tensor ¯τFictitious Kirchhoff stress tensor ΥjDissipative potentials φ∗Push-forward operation φ∗Pull-back operation Ψ Helmholtz free-energy function Ψvol Volumetric part of the strain energy function Ψiso Isochoric part of the strain energy function Ψ∞ vol Strain energy function describing the volumetric elastic response as t→ ∞ Ψ∞ iso Strain energy function describing the isochoric elastic response as t→ ∞ Symbol Description Ψiso jIsochoric part of the strain energy function responsible for the jrelaxation process Ω0Initial or reference configuration Ω Current or deformed configuration a0,g0Fibre directions in the reference configuration a,gFibre directions in the deformed configuration bLeft Cauchy-Green tensor ¯ bModified left Cauchy-Green tensor CRight Cauchy-Green tensor ¯ CModified right Cauchy-Green tensor C Elasticity tensor in material description C vol Volumetric part of the elasticity tensor in material description C iso Isochoric part of the elasticity tensor in material description c Elasticity tensor in spatial description c vol Volumetric part of the elasticity tensor in spatial description c iso Isochoric part of the elasticity tensor in spatial description ¯ c ∗Fictitious elasticity tensor in spatial description Dint Internal dissipation DCompressibility constant f Fat volume ratio F(t) Force FDeformation gradient tensor ¯ FModified deformation gradient tensor G(t) Relaxation function ¯ G(t) Reduced relaxation function giProny constants IaInvariants of C, with a= 1,2,3 IaPseudo-invariants of Cand a0⊗a0or g0⊗g0, with a= 4,5,6,7,8 ¯ IaModified invariants of C, with a= 1,2,3 ISecond-order unit tensor I Fourth-order unit tensor JVolume ratio J(t) Creep function pHydrostatic pressure ˜pModified hydrostatic pressure PFirst Piola-Kirchhoff stress tensor P Projection tensor QRotation tensor QjNon-equilibrium stress SSecond Piola-Kirchhoff stress tensor ¯ SFictitious second Piola-Kirchhoff stress tensor Svol Volumetric part of the second Piola-Kirchhoff stress tensor Siso Isochoric part of the second Piola-Kirchhoff stress tensor S∞ vol Volumetric part of the second Piola-Kirchhoff stress tensor when t→ ∞ S∞ iso Isochoric part of the second Piola-Kirchhoff stress tensor when t→ ∞ Siso jIsochoric part of the second Piola-Kirchhoff stress tensor corresponding to Ψiso j ¯ SjFictitious second Piola-Kirchhoff stress tensor corresponding to the jrelaxation process Symbol Description tTime TeElastic response function UDisplacement field vector in the material description uDisplacement field vector in the spatial description Contents 1 Introduction 11 1.1 Whystudyingthebreast?...................................... 11 1.2 Motivation .............................................. 15 1.3 Objectiveofthethesis........................................ 15 1.4 Contentsofthework......................................... 16 2 The female breast 19 2.1 Aglobalpicture............................................ 19 2.2 Adiposetissue............................................. 21 2.3 Glandulartissue ........................................... 23 2.4 Connectivetissue........................................... 24 2.5 Skin .................................................. 24 2.6 Methods to determine the fat proportion in the breast . . . . . . . . . . . . . . . . . . . . . . 25 2.7 Breastcancer............................................. 26 3 Basic concepts of Nonlinear Continuum Mechanics 29 3.1 Introduction.............................................. 29 3.2 Hyperelasticity ............................................ 30 3.2.1 Introduction ......................................... 30 3.2.2 Isotropic hyperelastic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2.3 Incompressible hyperelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.4 Compressible hyperelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2.5 Hyperelasticmodels ..................................... 34 3.2.6 Elasticitytensors....................................... 43 3.2.7 Fibre-reinforced materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3 Viscoelasticity............................................. 48 3.3.1 Quasi-linear viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3.2 Viscoelastic model with internal variables . . . . . . . . . . . . . . . . . . . . . . . . . 50 4 Mechanical properties of breast tissues 53 4.1 Introduction.............................................. 53 4.1.1 Experimental estimation of properties . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.1.2 Computational estimation of properties . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1.3 Mechanical properties of breast tissues. State of the art . . . . . . . . . . . . . . . . . 57 4.1.4 Mechanical properties of skin. State of the art . . . . . . . . . . . . . . . . . . . . . . 58 4.1.5 Mechanical properties of muscle. State of the art . . . . . . . . . . . . . . . . . . . . . 59 4.1.6 Cooper’sligaments...................................... 59 4.2 Model for the mixture of fat and glandular tissue . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2.1 StudyA............................................ 61 4.2.2 StudyB............................................ 65 4.2.3 StudyC............................................ 66 4.2.4 StudyD............................................ 68 4.2.5 Some remarks on the homogenized model . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3 Experimental characterisation of breast and abdominal tissues . . . . . . . . . . . . . . . . . 70 4.3.1 Testprotocol ......................................... 70 4.3.2 Data fitting algorithm: QLV model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.3 Data fitting algorithm: internal variables model . . . . . . . . . . . . . . . . . . . . . . 78 4.3.4 Description of the experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.3.5 Results ............................................ 82 4.3.6 Comparison between the Samani’s constants and the ones obtained here . . . . . . . . 98 5 Finite element models of the breast 101 5.1 Introduction..............................................101 5.2 FEmodelconstruction........................................103 5.3 Superficial3Dscanimages......................................109 5.4 Boundaryconditions.........................................110 5.4.1 Adiposetissue ........................................110 5.4.2 Muscle.............................................111 5.4.3 Glandulartissue .......................................112 5.4.4 Skin ..............................................112 5.5 Materialmodelsused.........................................112 5.6 FEModels ..............................................113 5.6.1 Muscle fixed to the ribs and sternum, without skin. . . . . . . . . . . . . . . . . . . . . 113 5.6.2 Muscle fixed to the ribs and sternum, with skin . . . . . . . . . . . . . . . . . . . . . . 119 5.6.3 All the boundaries of the muscle fixed, with skin . . . . . . . . . . . . . . . . . . . . . 119 5.6.4 Sensitivity of the model to the stiffness of the fat-gland properties, with skin . . . . . 121 5.6.5 Sensitivity of the model to the stiffness of the fat-gland properties, without skin . . . 123 5.6.6 Comparison between the material models based on Samani’s constants and the one obtainedinthisthesis ....................................125 5.6.7 Tissue layer in the interface between the muscle and breast tissues . . . . . . . . . . . 126 5.6.8 Complete lack of connection between the muscle and the breast tissues . . . . . . . . . 128 5.6.9 Spring connection between the muscle and the breast tissues . . . . . . . . . . . . . . 130 5.6.10 Analysis of the influence of the skin properties . . . . . . . . . . . . . . . . . . . . . . 131 5.6.11 Contact surface between the muscle and the breast tissues . . . . . . . . . . . . . . . . 136 5.6.12 Model with anisotropic properties for the breast tissues . . . . . . . . . . . . . . . . . 143 5.6.13Othermodels.........................................146 6 Closure 157 6.1 Summaryandconclusions......................................157 6.2 Originalcontributions ........................................159 6.3 Limitations ..............................................160 6.3.1 Experimentalwork......................................160 6.3.2 Computationalwork.....................................160 6.4 Futurework..............................................161 A Integrals for the QLV fitting algorithm 163 A.1 Polynomial strain energy function with five terms . . . . . . . . . . . . . . . . . . . . . . . . . 165 A.2 First order Ogden strain energy function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 A.3 Isotropic GOH strain energy function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 A.4 Combination of a neo-Hookean strain energy function and an exponential one . . . . . . . . . 171 B Safety manual 173 Bibliography 183 16 Introduction us diseased breasts or simulation of surgical procedures. In order to develop a model of the breast for any application to help clinicians in their daily work, it is crucial to understand how the breast behaves internally in the normal conditions, to be able to extend the model to more complex cases. Once the healthy breast FE model is constructed and validated, it is possible to add characteristics or change the properties of the tissue, to see how this new model behaves and to compare it to the healthy breast model, for instance, a breast with a tumour inclusion. It can also be helpful in the field of image analysis to track certain points in the images or to compare images in different positions or conditions. Regarding the breast reconstruction, having a validated healthy breast model is the first step in the process of simulating the surgery. Subsequently, the healthy model can be modified to reproduce the reconstructive surgery and it is also possible to optimize some parameters in the model (quantity of adipose tissue, shape, skin tension, etc.) in order to make the behaviour similar to the healthy one. To develop a FE model of the healthy breast, it is key to correctly characterise the mechanical properties of the breast tissues, which is then another important objective of this work. The abdominal adipose tissue will also be characterised, what can be useful for the reconstructive surgery simulation and, in general, for any other simulation conducted in the abdomen. Moreover, assuming the properties assessed here for the adipose tissue of the abdomen as valid for other regions, the applications could be quite numerous. 1.4 Contents of the work To develop the FE model of the healthy breast it is necessary to define its geometry, the load and displacement boundary conditions, the mechanical properties of the tissues involved and the interaction between these tissues. Therefore, this thesis can be divided in three big parts: an anatomical study of the structures and tissues which compose the breast, an experimental work in which the properties of the tissues are determined and a computational work in which the FE model is constructed and the simulations carried out. All these parts are explained subsequently divided into chapters. In chapter 2, the anatomy of the female breast is shortly presented. The different types of tissue which compose it are described, as well as their distribution and structure. Some methods to determine the fat proportion in the breast are also included. Finally, breast cancer is shortly introduced. In chapter 3, some fundamental concepts of Nonlinear Continuum Mechanics needed to understand this work are presented. This chapter is not intended to be a comprehensive text. Only the most important concepts are reminded. Some mathematical developments which are used in the subsequent chapters are also presented. Chapter 4 is about the mechanical properties of breast tissues. First, the previous experimental studies are described. Then, a homogenization procedure for the mixture of fat and glandular tissue is presented. Finally, the experimental test conducted is explained, including the test protocol, the algorithm developed Contents of the work 17 to fit the experimental curves and the results obtained. Chapter 5 is devoted to the FE models of the breast. As in chapter 4, the most relevant literature is introduced first. Then, the construction of the model is explained as well as the method to validate it. Finally, the different models are presented, explaining the differences between them and their performance to obtain the real breast deformation. In chapter 6, the thesis is summarised and the main conclusions extracted. The original contributions of this thesis and the future works are also presented. Appendix A includes the procedure to obtain the analytical stress in the QLV formulation. It implies a high number of integrals which have not been included in the main text for clarity. Finally, in appendix B, the safety manual developed for the Biomechanics Laboratory of the Department of Mechanical Engineering is attached. This manual was mainly made to work with human tissue specimens. Chapter 2 The female breast In this chapter, the anatomy of the female breast is shortly presented. Knowing the inner structure of the breast, the tissues which compose it and their distribution is of great importance to understand how the breast is deformed. Moreover, this knowledge is the starting point to develop a computational model of the breast. When modelling a biological system in a computational framework, the modeller should make some assumptions and simplifications through the process, relating the boundary conditions, the mechanical properties of the materials, the connections between them, etc. It is very important to always keep in mind the real anatomy, in order to know at any point if the subsequent simplifications and assumptions are correct, or at least, suitable. Finally, breast cancer is shortly introduced. 2.1 A global picture Breasts or mammary glands are the distinctive characteristic of mammals. In humans, these organs cover the female chest, but are not normally developed in men. The female breast is an organ which is present in the women’s body with a clear purpose: the breastfeeding of their children. The breast is located overlaying the women’s chest. It extends inferiorly from the second or third rib to the inframammary fold, which is approximately about the level of the sixth or seventh rib. Laterally, it extends form the sternum to the anterior axillary line. The extension into the axilla, known as the glandular Tail of Spence, is highly variable between subjects. The posterior face of the mammary gland is delimited by the fasciae of the pectoralis major, serratus anterior, external abdominal oblique and rectus abdominis muscles. The anterior part is covered by the skin [25, 26, 27, 28]. The female breast is mainly composed of glandular tissue, fat, connective tissue and skin. There also exist nervous tissue innerving the breast and blood vessels irrigating it. The glandular tissue is formed by lobules which are connected to the nipple through the lactiferous ducts in a radial arrangement. These lobules are separated by connective tissue and some fat is intermingled between them. The rest of the fat 20 The female breast Figure 2.1: Global view of the female breast (adapted from [26]). is surrounding the glandular tissue. There is some more connective tissue in the interface between the fat and the muscles, between the fat and the skin, and as ligaments inside the breast which are called Cooper’s ligaments (suspensory ligaments). Finally, the skin is externally covering the organ. The mammary development is initially identical in males and females. It starts approximately in the 4th - 6th week of gestation with the appearance of the milk lines (ectodermal thickening) on the ventral area. It extends from the axilla to the medial thigh (figure 2.2). In humans, the proximal and distal parts of the milk lines become atrophied while the pectoral part is the one responsible of the genesis of the breast. Inherited disorders may appear in 1% of women, as the development of an additional breast (polymastia) or a supernumerary nipple (polythelia). They can also appear in men. Once the milk lines have been created, during the remainder of gestation, they start to proliferate forming the lactiferous ducts. The areola is created during the 5th month. At birth, the breast, identical in men and women, is composed of radially arranged mammary lobes connected to the nipple through lactiferous ducts. This nipple is a little hollow formed by some Montgomery glands. Just after birth, the breasts may be active (with secretion of milk) due to the mother’s hormones. Then the mammary activity is stopped until puberty. Shortly after birth, the nipple is everted due to the proliferation of the mesoderm, and the pigmentation of the areola increases. In men, after birth, there is little mammary development. The presence of testosterone leads to the Adipose tissue 21 Figure 2.2: Milk lines (adapted from [26]). involution of the mammary gland. However, in female, the breast continues developing due to the action of hormones. Three phases can be differentiated. In the first one, the stroma and glandular tissue of the breast slowly grow from the birth to the puberty, but no lobular growing takes place. Then, from puberty to approximately 20 years old, there exists a fast growth because of the development of fat and connective tissue, but also due to the thickening of the ductal system. This phase takes place under the effect of estrogens, growth hormone and prolactin. Finally, when pregnancy occurs the breasts reach their functional state and, under the influence of progesterone, estrogens, prolactin and lactogen, the development is completed. They become secretory alveoli and start to accumulate the necessary substances to generate milk in the postpartum period. When nursing is abandoned, the alveoli return to their non-functioning state. By the age of 40 years the mammary glands start to atrophy. Then the glandular lobes begin to be replaced by connective tissue and fat. 2.2 Adipose tissue Adipose tissue is a loose connective tissue in which the adipocytes are the main cellular components. Adipocytes play a fundamental role in energy homeostasis. To survive, the human body should assure a constant energy contribution, despite of the fact that the availability of nutrients from the exterior may be variable. To compensate the energy demand of the body when there is not enough energy intake, the adipose tissue is able to store energy in an efficient way. The energy is stored in the lipid droplets as triglycerides (triacylglycerols). When food intake is higher than waste (in terms of energy), the excess is stored in the adipocytes. The triglycerides are the most concentrated way of storing energy in the human body. Moreover, the adipocytes have another functions besides the storing one. They secrete paracrine and endocrine 22 The female breast Figure 2.3: Development of an adipocyte (adapted from [27]). substances to regulate the energy metabolism. The fat also has a structural function, serving as a padding to protect other organs. There are two types of adipose tissue: white (unilocular) and brown (multilocular). The names come from the aspect that the cells have under the microscope and the color of the tissue when it is fresh. The white adipose tissue is predominant in the human body. The brown one is present during the fetal stage, but its quantity decreases along the first ten years after birth. The unilocular adipose tissue forms a subcutaneous layer. It has an important isolating function. In the body it is normal to find big amounts of this adipose tissue in the abdominal region, the gluteal area, the axilla and the thigh. In both sexes the pectoral region is a preferential area to store fat. In fact, the fat is the main component of the non lactating breast. The origin of these adipocytes is an undifferentiated mesenchymal stem cell (figure 2.3). They have a spheroidal shape when they are isolated, but become an oval or a polyhedron when grouped. The big size of these cells is due to the storage of lipids, which causes the movement of the nucleus towards the border and its flattening. They are surrounded by reticular fibres (type III collagen) generated by them (figure 2.4). Moreover, the adipose tissue has a high blood supply, with many capillaries within its structure. The amount of white adipose tissue in the body is regulated by two physiological systems. The first one, related with the regulation of weight in the short-term, controls the appetite and the daily metabolism. The second one, related with the regulation of weight in the long-term, controls the appetite in a continuous way Glandular tissue 23 Figure 2.4: Adipocytes (adapted from [27]). (during months or years). The cells of the multilocular adipose tissue are smaller than those of the unilocular one. The nucleus of the multilocular adipocyte is eccentric but not flattened as in the unilocular one. They also derive from undifferentiated mesenchymal stem cells (figure 2.3). In the breast, the fat (as said before, white adipose tissue) is surrounding the glandular tissue and filling the spaces between the lobes. The amount of fat in the breast is highly variable. It is responsible for the main part of the volume, shape and consistency of the breast. The percentage of body fat has an influence on the amount of adipose tissue in the breast, although it seems that some women, independently of their body fat, are more prone to accumulate fat in their breasts. The breasts with a higher percentage of fat are softer, while those with a higher percentage of glandular tissue are normally firmer and nodular. The amount of fat in the breast is increased once the glands start to atrophy after lactation and menopause. 2.3 Glandular tissue Mammary glands are modified apocrine sweat glands (sudoriparous glands) which are developed due to the action of sex hormones. The mature inactive breast is composed of 15 −20 irregular lobes of ramified tubuloalveolar glands. The lobes, separated from each other by fibrous connective tissue, have a radial distribution starting from the nipple. Around 20 main lactiferous ducts (2 −4.5mm in diameter) drain the mammary lobes to the nipple. Each lobe is subdivided into little lobes (200 −600 µm ), called terminal duct lobular units (TDLU) or lobules, composed of hundreds of secretory acini each one of which is connected to the lactiferous ducts through individual interlobular ducts (40 −60 µm in diameter) [27, 29, 30]. All the system is inside the connective stroma and the adipose tissue. The nipple and the areola are highly pigmented in contrast to the surrounding skin and a little bit wrinkled. Beneath the skin, there are smooth muscle fibre bundles arranged radially and circumferentially, and in a longitudinal direction along the lactiferous ducts. These muscle fibres are responsible for the erection 24 The female breast of the nipple in response to different stimuli. The areola has sebaceous glands, sweat glands and modified mammary glands (Montgomery glands). These glands, which have an intermediate structure between the sweat glands and the mammary glands, produce little lumps in the skin of the areola. The lactiferous ducts end in the nipple through a small pore. Under the areola, each duct has a dilated portion called lactiferous sinus. When the mammary glands are inactive, the glandular component is limited and it mainly consists of ducts. During menstruation, and principally when the ovulation is about to start, secretory cells get bigger and ducts get wider, due to the action of oestrogen. During pregnancy, glandular tissue suffers some changes and fat and connective tissue are decreased. 2.4 Connective tissue In the breast there are several structures made from connective tissue. The breast is developed and contained inside some supporting layers of the superficial fascia. The superficial layer of the superficial fascia is located near the skin, and is not always differentiated from it. The superficial layer of the superficial fascia is the external layer which covers the breast parenchyma. A variable amount of fat can be found between the dermis and the parenchyma. The deep layer of the superficial fascia is on the posterior surface of the breast. Some loose connective tissue is between the deep layer of the superficial fascia and the deep fascia which is covering the pectoralis mayor muscle, the rectus abdominis muscle, the serratus anterior muscle and the external abdominal oblique muscle (see figure 2.1). The deep layer of the superficial fascia is pierced by some fibrous supports, called Cooper’s ligaments, and by blood vessels, nerves and lymph vessels. Cooper’s ligaments connect the deep layer of the superficial fascia and the skin which covers the breast, passing through all the breast. These suspensory ligaments are connected to the deep muscle fascia. Due to the fact that these ligaments are not tense, they allow the movement of the breast over the deep muscle fascia. These ligaments can be stretched and elongated because of pregnancy, ageing or fluctuations in weight. The loss of the stiffness of these connective tissues contributes to the ptosis of the breast. 2.5 Skin The skin is the biggest organ in the body. It is composed of two main layers (figure 2.5): •The epidermis, formed by a keratinized stratified epithelium which is constantly growing though it always maintains its normal thickness due to the desquamation process (skin peeling). •The dermis, composed of a dense connective tissue to provide mechanical support. Methods to determine the fat proportion in the breast 25 Figure 2.5: Layers of the skin. The skin is involved in many functions and processes of the human body: •It is a protection barrier against physical, chemical or biological agents from the surroundings. •It provides immune information. •It takes part in the body temperature homeostasis by regulating the heat exchange and the loss of water. •It sends information about the surroundings obtained through the senses. •It has endocrine functions. •It plays a role in the excretion. The thickness of the skin varies from less than 1mm to more than 5mm depending on the body region. Obviously, there are two regions in the body where the skin is clearly different: the palm of the hands and the sole of the feet. These regions are subjected to intense friction, they do not have hair and they have a much thicker epidermis than in any other part of the body. In the breast, the quality, thickness and elasticity of the skin have a big influence in its appearance. The skin is influenced by hormonal changes, weight variations, gravity and ageing, which can cause a reduction in the thickness or a loss of elasticity. Stretch marks normally appear after a weight loss or a labour. 2.6 Methods to determine the fat proportion in the breast Knowing the fat proportion in the breast is important for many reasons. The main one is based on the fact that breast tissue density has been identified as an important risk factor for cancer development and its analysis is useful for quantitative mammographic assessment. Moreover, for many surgeries it is useful to have an estimation of the fat proportion beforehand. This is linked to the availability of assessing the fat 32 Basic concepts of Nonlinear Continuum Mechanics P= 2F∂Ψ(C) ∂C,(3.7a) S=2∂Ψ(C) ∂C= 2 ∂Ψ ∂I1 ∂I1 ∂C+∂Ψ ∂I2 ∂I2 ∂C+∂Ψ ∂I3 ∂I3 ∂C= 2 ∂Ψ ∂I1 +I1 ∂Ψ ∂I2I−∂Ψ ∂I2 C+I3 ∂Ψ ∂I3 C−1.(3.7b) Recalling the relation between the Cauchy stress tensor and the second Piola-Kirchhoff stress tensor, σ=J−1FSFT, with J= det(F) the volume ratio, and using equation (3.7b) and the definition of the left Cauchy-Green tensor b=FFT, it is possible to express the Cauchy stress tensor in terms of the three invariants: σ= 2J−1I3 ∂Ψ ∂I3 I+∂Ψ ∂I1 +I1 ∂Ψ ∂I2b−∂Ψ ∂I2 b2.(3.8) 3.2.3 Incompressible hyperelasticity In practice, there are many processes in which the materials involved behave as incompressible or at least very close to it (nearly incompressible or quasi-incompressible). In the case of biological tissues, as many of them are mainly composed of water, and water is an incompressible material, they can be treated as incompressible as well. Modelling incompressible materials can lead to numerical problems that may compromise the convergence of the problem. This is the reason why it is advisable to model the incompressible materials as quasi-incompressible, that is to say, allowing a certain volumetric deformation in order to facilitate convergence. The restriction which characterises incompressibility is J= 1. A material subjected to an internal constraint is called a constrained material. In this case, the strain energy function can be written as: Ψ = Ψ(F)−p(J−1) ,(3.9) where pcan be seen as a Lagrange multiplier, which can be identified with the hydrostatic pressure. The scalar pcan only be determined from the equilibrium equations and the boundary conditions. Equation (3.7b) and (3.8) can be rewritten for an incompressible isotropic hyperelastic material as: S= 2∂Ψ(I1, I2) ∂C−∂[p(I1/2 3−1)] ∂C=−pC−1+ 2 ∂Ψ ∂I1 +I1 ∂Ψ ∂I2I−2∂Ψ ∂I2 C,(3.10) σ=−pI+ 2 ∂Ψ ∂I1 +I1 ∂Ψ ∂I2b−2∂Ψ ∂I2 b2.(3.11) Hyperelasticity 33 3.2.4 Compressible hyperelasticity A material which can experience volumen changes is said to be a compressible material. In numerical simulations, quasi-incompressible materials are usually treated with the same formulation as compressible ones to ensure an easier and faster convergence. Most materials behave in a different way under bulk and shear deformations, so that it is convenient to split the strain energy function in the so-called volumetric and isochoric parts. In particular, Fand Care decomposed into volume changing (dilatational) and volume preserving (distortional) parts: F= (J1/3I)¯ F=J1/3¯ F,C= (J2/3I)¯ C=J2/3¯ C,(3.12) where J1/3Iand J2/3Iare associated with the change of volume, and ¯ Fand ¯ C=¯ FT¯ F, called the modified deformation gradient tensor and modified right Cauchy-green tensor respectively, represent the distortion deformation. Then, it is possible to express the strain energy function with a decoupled representation: Ψ(C)=Ψvol(J)+Ψiso(¯ C).(3.13) The multiplicative decomposition of tensor Fis supported in the field of computational mechanics if working with quasi-incompressible (or slightly incompressible) materials to avoid numerical complications by treating separately the volumetric and isochoric parts. In case of using a standard displacement based method (Galerkin method) in which only the displacement field is discretized, numerical problems can arise if modelling nearly incompressible materials due to locking phenomena (in general, ill-conditioning of the stiffness matrix with over-stiffening of the system). Because of this, multi-field variational principles and a mixed finite element method (FEM) are used. In these methods, not only the displacement field is discretized, but also other fields like, for instance, the internal pressure or the volume ratio. The idea is to discretize variables in an independent way in order to achieve nonlocking and stable numerical solutions. The second Piola-Kirchhoff stress tensor Scan also be split from equation (3.7b): S= 2∂Ψvol(J) ∂C+ 2∂Ψiso(¯ C) ∂C=Svol +Siso ,(3.14a) Svol =JpC−1,Siso =J−2/3 I −1 3C−1⊗C:¯ S,(3.14b) where I is the fourth-order unit tensor. The hydrostatic pressure pis defined as: p=dΨvol(J) dJ(3.15) 34 Basic concepts of Nonlinear Continuum Mechanics and the tensor ¯ S, called the fictitious second Piola-Kirchhoff stress tensor, as: ¯ S= 2∂Ψiso(¯ C) ∂¯ C.(3.16) It is important to remark that, unlike incompressible materials, in which the hydrostatic pressure is given by equilibrium equations and boundary conditions, in compressible materials the scalar function pis defined by a constitutive equation. The Cauchy stress tensor is also split into the volumetric and isochoric parts, σvol yσiso, in the following way: σvol = 2J−1b∂Ψvol(J) ∂b=pI,(3.17a) σiso = 2J−1b∂Ψiso(¯ b) ∂b= I −1 3I⊗I:¯σ= P :¯σ= dev¯σ,(3.17b) where P = I −1 3I⊗Iis the projection tensor, a fourth order tensor which performs the transformation of a second order tensor into its deviatoric part. The fictitious Cauchy stress tensor, ¯σ, is defined as: ¯σ= 2J−1∂Ψiso(¯ b) ∂¯ b¯ b,(3.18) with ¯ b=¯ F¯ FTthe modified left Cauchy-Green tensor. Finally, the decoupled strain energy function can be expressed as a function of the invariants of ¯ Cor ¯ b, for compressible isotropic hyperelastic materials, similarly to equation (3.6): Ψ = Ψvol(J)+Ψiso[¯ I1(¯ C),¯ I2(¯ C)] = Ψvol(J)+Ψiso[¯ I1(¯ b),¯ I2(¯ b)] .(3.19) 3.2.5 Hyperelastic models Many hyperelastic models have been proposed in the literature, both for rubber-like materials and for biological materials. These models arise from the need to reproduce the behaviour of real materials. Normally, to characterise the mechanical behaviour of a given material, a set of experimental tests is carried out. Once the results are obtained, a model is needed to fit the results. The phenomenological approach proposes strain energy functions which try to fit the data in a purely mathematical way, that is to say, they do not take into account the internal structure of the material, they just treat it as a continuum. As a consequence, the material parameters which are present in these models do not usually have a physical meaning. This is a fast and easy procedure, but it loses touch with reality and physical sense. On the other hand, other Hyperelasticity 35 approaches try to model the inner structure to some extent, giving some physical meaning to the material parameters. The good point is that, theoretically, it is possible to obtain some of the parameters from the study of the material structure, being closer to reality and reducing the number of parameters which need to be fitted in a mathematical way. The problem is that normally it is a hard task to obtain such a model which really works. Moreover, only one type of experimental test is generally fitted, for instance, a uniaxial tension test. This means that the computational model will behave as expected for this loading case, but not necessarily well for another type of load, shear stress for example. Many models fail to reproduce the behaviour of a material subjected to several loading conditions. Thus, it is important to know which material model is being used, how it behaves subjected to different loading conditions and which is the stress state the body is going to suffer in practice. In this section, some strain energy functions will be presented, and in particular, those that were used in this work. 3.2.5.1 Incompressible models Polynomial models. Probably the simplest model in the literature is the neo-Hookean model. The strain energy function of this material is: Ψ = C10(I1−3) ,(3.20) where C10 is a stress-like parameter. The simplicity of this model resides in its dependence on only one material parameter and the first invariant I1. Due to this fact, it has been extensively used in the literature, for example in the simulation of the mechanical behaviour of the breast to predict deformations [18] or to simulate the behaviour of the ligaments of the knee [37]. This strain energy function is also used to model the behaviour of the isotropic matrix in the well known Holzapfel-Gasser-Ogden (HGO) fibre-reinforced model [38]. Another simple strain energy function, proposed at first by Mooney and reformulated by Rivlin, is the Mooney-Rivlin model: Ψ = C10(I1−3) + C01(I2−3) .(3.21) Chen et al. [39] used this model to simulate the behaviour of the articular disc of the temporomandibular joint. In general, both of them belong to the family of the polynomial models. The polynomial strain energy functions have the following general form: 36 Basic concepts of Nonlinear Continuum Mechanics Ψ = N X i+j=1 Cij(I1−3)i(I2−3)j,(3.22) where Cij are stress-like parameters. Thus, from equation (3.22), many different strain energy functions can be obtained varying the number of terms. In particular, the neo-Hookean and the Mooney-Rivlin models presented in equations (3.20) and (3.21) respectively, are particular cases of this polynomial function. In the case of the neo-Hookean model for only one term, and in the case of the Mooney-Rivlin model for N= 1. For example, Li et al. [40] used a 3 terms polynomial strain energy function to simulate the fibrocartilage of the human pubic symphysis. Samani and Plewes [41] modelled the properties of the adipose tissue and the glandular tissue of the breast with a 5 terms polynomial strain energy function (N= 2). Ψ = C10(I1−3) + C01(I2−3) + C11(I1−3)(I2−3) + C20(I1−3)2+C02(I2−3)2.(3.23) Let us consider an isotropic incompressible hyperelastic material, whose behaviour is modelled by means of a polynomial function with N= 2, subjected to uniaxial tension in the zaxis. This load case will be of interest in following sections, and that is the reason why the theoretical development is presented here. The Cauchy stress tensor in this case can be written: σ=       0 0 0 0 0 0 0 0 σT z an       .(3.24a) The aim is to obtain the analytical expression of σT z an as a function of the deformation. For uniaxial compression, the process is equivalent just by substituting σC z an for σT z an. The subscript “an” stands for analytical solution. The deformation gradient tensor is given as a function of the principal stretches λa: F=       λ10 0 0λ20 0 0 λ3       .(3.24b) Using the incompressibility constraint, J= 1, and the symmetry of loads and displacement boundary conditions, which imply that the stretches in both axes xand ymust be equal, λ1=λ2, it is possible to obtain the relation between the principal stretches: Hyperelasticity 37 J= 1 = det(F) = λ1λ2λ3,(3.24c) If λ3is termed simply λ, equation (3.24c) implies: λ1=λ2=1 √λ.(3.24d) Therefore, the deformation gradient tensor in (3.24b) has the form: F=       1 √λ0 0 01 √λ0 0 0 λ       (3.24e) and the left Cauchy Green tensor b: b=       1 λ0 0 01 λ0 0 0 λ2       .(3.24f) Using equation (3.11) the boundary conditions yield: σT y an =−p+ 2hC10 +C01I1+C11I1I1−3+C11I2−3+ 2C20I1−3+ + 2C02I1I2−3i1 λ−2hC01 +C11I1−3+ 2C02I2−3i1 λ2= 0 ,(3.25a) from which the hydrostatic pressure pcan be worked out: p=2 λ42C02(−1 + λ)2(1 + 2λ+λ3+ 2λ4) + λλC01 + 2C20 +C10λ− −3C20λ+ (C01 +C20)λ3+C11(−1 + λ)2(3 + 3λ+ 2λ3+λ4).(3.25b) 38 Basic concepts of Nonlinear Continuum Mechanics 0.5 1 1.5 −40 −30 −20 −10 0 10 20 30 λ Cauchy stress (kPa) 1 term polynomial function 3 terms polynomial function 5 terms polynomial function Figure 3.1: Stress-stretch curve for uniaxial tension-compression for the 1, 3 and 5 terms polynomial models with the parameters: C10 = 0.3 kPa, C01 = 0.3 kPa, C11 = 2.25 kPa, C20 = 3.8 kPa and C02 = 4.72 kPa. These constants were proposed by [41] for the breast adipose tissue. Finally, the Cauchy stress tensor in zdirection is calculated after substituting p: σ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.(3.26) The Cauchy stress for the neo-Hookean and Mooney-Rivlin materials, and in general, for any polynomial function with a number of terms less than 5, can be derived from equation (3.26) by cancelling the corresponding terms. For example, in the neo-Hookean model C01 =C11 =C20 =C02 = 0 and in the Mooney-Rivlin model C11 =C20 =C02 = 0. In figure 3.1 an example of the Cauchy stress-stretch (σ−λ) curve in a uniaxial test is shown for the polynomial models with one (neo-Hookean), three and five terms. It can be noticed the different behaviour of the three models and how the curve becomes more markedly non-linear as the number of terms of the polynomial strain energy function increases, in compression (λ < 1) as well as in tension (λ > 1). In fact, for the stretch range presented in the graph and for the constants used, it is possible to notice that the behaviour of the neo-Hookean material is very close to linear. As more terms are added (3 and 5 terms polynomial functions) the curves become stiffer, more clearly in compression. Another load case, also used later in this work, is a material subjected to a pure shear stress state. Let Hyperelasticity 39 Figure 3.2: Pure shear stress state. us consider an isotropic incompressible hyperelastic material, with its mechanical behaviour described by a polynomial strain energy function with N= 2 and with a stress state shown in (3.27a). σ=       0 0 0 0 0 σyz 0σzy 0       .(3.27a) The aim, as in the previous case, is to obtain the stresses as a function of deformation. However, implementing this stress state in a FE code causes some problems when the displacement boundary conditions are applied. The classical way of representing a pure shear stress state in linear elasticity theory is shown in figure 3.2. If the displacement boundary conditions are applied in this way in a large displacements and strains problem, other components of the stress tensor apart from σyz and σzy will appear. This problem is solved using the fact that this stress state can be expressed in the principal directions of the stress tensor (at 45ofrom the yand zdirections shown in figure 3.2) as a tension and a compression of the same magnitud, that is: σppl =QTσQ =       0 0 0 0σS y an 0 0 0 σS z an       ,(3.27b) where Qis a rotation matrix. Now, the displacement boundary conditions can be easily applied. In the linear elastic case, the behaviour of the material in compression and tension is equivalent, so that if σyz =σzy =σ, it turns out that σS y an =−σS z an =σ. In the following mathematical development, applied to a hyperelastic material subjected to large strains and displacements, the behaviour in compression and tension may be different, in contrast to the linear elastic case. Therefore it is not true in general that σS y an =−σS z an. Instead, the stress tensor is characterised by two different components. The aim is to obtain the stresses 40 Basic concepts of Nonlinear Continuum Mechanics as a function of deformation. To do that, the calculation of the deformation gradient is the starting point. Imposing zero displacements in xdirection (perpendicular to the yz plane in figure 3.2), and keeping in mind that this is an incompressible material (equation (3.24c)), the cross-sectional area of the probe in the yz plane must remain constant during deformation. It can be proved that accomplishing these conditions and applying the tension in zdirection and the compression in ydirection, the tensor Fis: F=       1 0 0 0 1/λ 0 0 0 λ       (3.27c) and the left Cauchy-Green tensor b: b=       λ20 0 01 λ20 0 0 λ2       .(3.27d) Using equation (3.11) the boundary yield: σS x an =−p+ 2C10 −6C20 +C01−1 + I1+2C20 +C11−4 + I1I1+ + 2C02−1 + I1−3 + I2) + C11I2= 0 ,(3.28a) from which the hydrostatic pressure pcan be worked out: p=2 λ4C11 + (C01 −C11 + 2C20)λ2+ (C10 −4C20)λ4+ + (C01 −C11 + 2C20)λ6+C11λ8+ 2C02(−1 + λ2)2(1 + λ4).(3.28b) Finally, the components of the Cauchy stress tensor, σS y an and σS z an can be expressed as: σS y an =−2(−1 + λ2) λ4C11 + 2C20 + (C10 −C11 −4C20 + 2C02)λ2+ + (C01 −C11 + 2C20 −4C02)λ4+ (C11 + 2C02)λ6,(3.29a) Hyperelasticity 41 σS z an =2(−1 + λ2) λ4C11 + 2C02 + (C01 −C11 + 2C20 −4C02)λ2+ + (C10 −C11 −4C20 + 2C02)λ4+ (C11 + 2C20)λ6.(3.29b) In figure 3.3, an example of both components of the Cauchy stress tensor is shown for the polynomial models with one (neo-Hookean), three and five terms. As in the previous uniaxial example, it can be observed that the more terms the polynomial function has, the stiffer the response is. Again, the expression for the neo-Hookean and Mooney-Rivlin materials, and in general for any other polynomial model with less terms can be derived from equations (3.29a) and (3.29b), by cancelling the corresponding parameters: C01 =C11 =C20 =C02 = 0 for the neo-Hookean and C11 =C20 =C02 = 0 for the Mooney-Rivlin. Rubin-Bodner model. The Rubin-Bodner model, later used in this work, was proposed by Rubin and Bodner for biological tissues [42]. Then it was used by Barbarino et al. [43] in a FE model of the face soft tissues. In particular, these authors used this strain energy function to model the skin and the muscles. They did not consider the transient and dissipative behaviour (viscoelastic), so that only the long term hyperelastic terms were included. The Rubin-Bodner strain energy function is: Ψ = µ0 2q(eqg −1) ,(3.30) with: g=g1+g2= 2 m1J−1−ln(J)+ (1 −ω)m2(¯ I1−3) ,(3.31) where ¯ I1is the modified first invariant, g1characterises the response of the fully elastic dilatation and g2 the fully elastic distortion, m1and m2are material parameters and ωis an aging term. 3.2.5.2 Compressible and quasi-incompressible models Analytically, the incompressibility condition is treated as explained in the previous sections. In the FEM, incompressibility can lead to several numerical problems due to the fact that the stiffness matrix becomes singular. To avoid these problems, materials are treated as quasi-incompressible, splitting the strain energy function in a volumetric and an isochoric part. But even using quasi-incompressibility, as explained in section 3.2.4, if a standard displacement based method is used, some numerical difficulties appear due to the fact that the stiffness matrix is almost singular because the effective bulk modulus is very big compared to the effective shear modulus. That depends on how strictly the condition is imposed, i.e., how close is 48 Basic concepts of Nonlinear Continuum Mechanics the abdominal adipose tissue [47]. The proposed strain energy function is: ¯ Ψ = C10 2(¯ I1−3) + k1 2k2X i=4,6ek2[κ¯ I1+(1−3κ)¯ Ii−1]2−1.(3.58) In the above equation, the parameters C10,k1and k2have the same meaning that in the original model. The new parameter κ∈[0,1/3] accounts for the fibre dispersion. If its value is 0, the fibres are perfectly oriented in one direction without dispersion, whereas if its value is 1/3, the fibres are uniformly distributed in all directions and the material becomes isotropic. 3.3 Viscoelasticity The majority of biological materials, and many other materials used in different applications, are characterised by constitutive equation which depends explicitly on time. If, besides, there is no damage, the behaviour is called viscoelasticity. The viscous part is responsible for the dependence on time (strain rate) and for the dissipation of energy during a loading cycle. Unlike elastic materials, viscoelastic ones exhibit hysteresis. The elastic part is responsible for the dependence on the strain and for finally recovering the initial state when all the loads are removed. Another feature of viscoelastic materials is the creep and stress relaxation phenomena. Creep is revealed when, after applying a constant stress, the material continues deforming beyond the onset of the stress. The stress relaxation is the opposite: the decrease in stress in a material subjected to a constant deformation. Theoretically, in a creep test (or in a stress relaxation test) the stress (deformation) is applied as a step function. In practice, it is impossible to apply a step-like load. Therefore, the step must be approximated by an initial ramp with a slope as steep as possible. In a creep test under a unit stress, the temporal evolution of strain is called creep function, J(t). In a stress relaxation test under a unit strain, the temporal evolution of stress is known as relaxation function, G(t). 3.3.1 Quasi-linear viscoelasticity Quasi-linear viscoelasticity was first introduced by Fung [48] as a simplified way to deal with viscoelasticity. It allows characterizing the material behaviour as time dependent (viscoelasticity) and using a nonlinear stress-strain relationship. It assumes that the stress produced by a deformation step function can be written as: σ(λ, t) = G(t)Te(λ),(3.59) where λis the stretch produced by the deformation and G(t) is the reduced relaxation function. This is a dimensionless function normalized such that G(0) = 1. The function Te(λ), called the elastic response function, is expressed in the same units as the stress and represents the instantaneous stress which appears Viscoelasticity 49 in the material at the onset of the stretch. Linearity is maintained as long as the Boltzmann’s superposition principle is applied to give the stress evolution as: σ(t) = Zt −∞ G(t−τ)dTe[λ(τ)] dτ dτ . (3.60) Non-linearity is introduced in the formulation through the function Te(λ). This mixture of linear and non-linear formulations is what gives the name to the quasi-linear viscoelastic model. Considering that the material has initially no residual stresses or deformations: σ(t) = Zt 0 G(t−τ)dTe[λ(τ)] dτ dτ ∀t∈[0,∞).(3.61) It is possible to apply the chain rule to simplify the above equation, for subsequent calculations, as will be seen in later chapters: σ(t) = Zt 0 G(t−τ)dTe(λ) dλ dλ(τ) dτ dτ ∀t∈[0,∞).(3.62) The elastic response function, Te(λ), can be derived from any strain energy function, like for example, the neo-Hookean or the Mooney-Rivlin, as the stress in a uniaxial test for a given stretch λ. The Prony series are commonly used for the relaxation function, G(t): G(t) = E∞+ n X i=1 Eie−t/τi,(3.63) with nrepresenting the Maxwell units in a generalised Maxwell model, Eiand τiare, respectively, the stiffness and the relaxation time of each unit. The reduced relaxation function, G(t), is obtained from G(t) by normalizing it such that G(0) = 1. Thus, it is rewritten in terms of the Prony constants, gi, as: G(t) = g∞+ n X i=1 gie−t/τi,(3.64) with gi=Ei E∞+ n P i=1 Ei and g∞+ n X i=1 gi= 1 .(3.65) Usually, the number of terms in the Prony series is fixed by the modeller a priori in order to obtain a good fit of the experimental data. In the fitting of experimental data with Prony series, it is important to take into account that different sets of τiand gican lead to equal goodness of fit, as shown by Lanczos [49]. 50 Basic concepts of Nonlinear Continuum Mechanics This means that the solution fitted may be not unique and obviously this is not convenient for statistical purposes. Troyer et al. [50] demonstrated that fixing the time constants τiensures that the mathematical form of equation (3.64) yields a unique solution, i.e, it is not possible to have two sets of githat fit the curve equally well. Because of this, it is a good practice to fix the values of τibeforehand and only fit the constants gi. 3.3.2 Viscoelastic model with internal variables The internal variables, often called hidden or history variables, are used to model materials which involve dissipation. Although these variables try to model a real characteristic of the material, they are not observable, therefore being treated in a purely phenomenological way. The thermodynamic state of a material depends on the deformation gradient, F, its temperature and the internal variables. In this thesis, the temperature is assumed constant. The equations governing this behaviour are just briefly presented. For a deeper insight the reader is referred to the text of Holzapfel [36]. In contrast with the previously presented viscoelastic model, which is only valid for a uniaxial stress state, the viscoelastic model presented in this section is designed for a 3D stress state. However, for a uniaxial stress state under small strains and displacements, both models coincide. The strain energy function, using the decoupling into the volumetric and isochoric parts, can be expressed as: Ψ=Ψ∞ vol(J)+Ψ∞ iso(¯ C) + m X j=1 Υj(¯ C,Γj),(3.66) where Ψ∞ vol and Ψ∞ iso are respectively the volumetric and isochoric elastic response when time tends to infinity. The last term in equation (3.66) represents the non-equilibrium state and is called the dissipative potential. Tensors Γjdetermine the relaxation/creep behaviour, with j= 1, ..., m the number of “viscoelastic processes” with their corresponding relaxation times τj, which play the same role as τiin the QLV model. The second Piola-Kirchhoff stress tensor takes the form: S=S∞ vol +S∞ iso + m X j=1 Qj,(3.67) with Qjrepresenting the non-equilibrium stresses, and S∞ vol and S∞ iso are respectively the fully elastic volumetric and isochoric contributions to the second Piola-Kirchhoff stress tensor. The evolution equations for these internal variables are taken, by analogy with the non-equilibrium stresses in a generalized Maxwell model, as: ˙ Qj+Qj τj =˙ Siso j,(3.68) Viscoelasticity 51 where Siso jis the isochoric second Piola-Kirchhoff stress tensor which corresponds to the strain energy function Ψiso j, responsible for the jrelaxation process and defined in the subsequent equations. It is defined with the same structure of the right term of equation (3.14b) and equation (3.16): Siso j=J−2/3 P :¯ Sj,¯ Sj= 2∂Ψiso j(¯ C) ∂¯ C,(3.69) being ¯ Sjthe fictitious second Piola-Kirchhoff stress. The solution of equations (3.68) for t ∈(0,T] is: Qj=e−T/τjQj0++Zt=T t=0+ e−(T−t)/τj˙ Siso j(t)dt , (3.70) where Qj0+is the instantaneous response at the initial conditions. If the following assumption is made: Ψiso j(¯ C) = β∞ jΨ∞ iso(¯ C),(3.71) the stress Siso jcan be simplified as: Siso j=β∞ jS∞ iso(¯ C).(3.72) The β∞ jconstants are dimensionless strain energy factors, which make possible the simplification of the model through equations (3.71) and (3.72). These β∞ jconstants are equivalent to the giconstants in the QLV model. Chapter 4 Mechanical properties of breast tissues 4.1 Introduction To develop a biomechanical FE model (and in general, any computational model) it is very important to characterise the behaviour of the materials which compose it, in order to obtain a reliable model. As seen in the previous chapter, the balance equations are valid for any material, but there are more unknowns than equations. So, it is necessary to add the constitutive equations of the material to solve the problem. In this chapter, the experimental and computational procedures to estimate the mechanical properties of the tissues are presented in this introduction, along with the state of the art of the mechanical characterization of the tissues which compose the breast. Then, in section 4.2, a model for the mechanical behaviour of the mixture of fat and gland is included. This study was published in the article “J.L. Calvo-Gallego, J. Mart´ınez-Reina, J. Dom´ınguez. A polynomial hyperelastic model for the mixture of fat and glandular tissue in female breast. International Journal for Numerical Methods in Biomedical Engineering 2015; 31(9). doi: 10.1002/cnm.2723”. Finally, the experimental tests carried out in this thesis to characterize the tissue are presented in section 4.3. A large part of the experimental procedure used in this section was published in the article “M.S. Commisso, J.L. Calvo-Gallego, J. Mayo, E. Tanaka, J. Mart´ınez-Reina. Quasi-Linear Viscoelastic Model of the Articular Disc of the Temporomandibular Joint. Experimental Mechanics 2016; 56:1169-1177”. Another part of the experimental procedure has been submitted as “J.L. Calvo-Gallego, M.S. Commisso, J. Dom´ınguez, E. Tanaka, J. Mart´ınez-Reina. Effect of freezing and storage time on the elastic and viscous properties of the disc of porcine temporomandibular joint” to the Journal of the Mechanical Behavior of Biomedical Materials. 54 Mechanical properties of breast tissues 4.1.1 Experimental estimation of properties The most convenient way to determine the properties of a biological material is to carry out experimental tests. The final objective is to obtain a set of parameters which characterise a given strain energy function. This function represents the behaviour of the material and provides a mathematical framework to work with. Thus, the normal process is as follows: first, the tissue specimens are extracted from the donor (human or animal) and they are tested following a certain protocol; once the results of the tests are obtained, a strain energy function should be chosen; then the non-linear Continuum Mechanics equations are used to determine an analytical solution for the stress-strain relation of the performed test, which will be fitted to the experimental results. The choice of the strain energy function will depend on the shape of the experimental curve and on the results of the fitting algorithm. Normally, it is necessary to use several different functions, checking their goodness of fit to choose the best. The output of the fitting algorithm is a set of material parameters which defines the behaviour of the material together with the strain energy function used. In the following, some parameters of the testing of living tissues and their influence on the results will be discussed. 4.1.1.1 Effect of freezing In experimental tests of biological tissues, there are many things that should be taken into consideration. For instance, the tissues should be tested submerged in saline solution at 37◦C, to maintain the physiological conditions of temperature and hydration, although some authors indicate that there is no difference between testing the specimens at 37◦C or at room temperature [51]. A great debate still exists regarding the influence of freezing before testing on the mechanical properties of biological tissues and the influence of the time elapsed from the excision of the specimens to the test. Many authors claim that the cold storage causes changes in the mechanical properties of the tissues, as for example Linde and Sørensen [52], who analysed the effect of different storage methods on the elastic and viscoelastic properties of human trabecular bone under unconfined compression cyclic loading. They did not notice changes in stiffness but a significant change in viscoelastic properties. Clavert et al. [53] performed uniaxial tensile relaxation tests on human tendons (biceps brachii tendons) comparing the response obtained with immediate testing with that obtained after freezing at −30◦C during 2 weeks, finding an altered ultimate tensile stress and Young’s modulus. Kennedy et al. [54] carried out cyclic indentation tests on bovine articular cartilage, comparing a group tested fresh with specimens frozen at −20◦C and −80◦C. They checked that freezing decreased the stiffness. Chow and Zhang [55], who performed biaxial tensile tests on bovine aortic tissue stored at different temperatures and during intervals of different lengths, found several changes in the mechanical properties due to cold storage. Ternifi et al. [56] studied the influence of storage temperature on porcine kidney with shear wave elastography, concluding that freezing causes a decrease in Introduction 55 the shear modulus of the superficial zone but not in the central one. Opposite to the previous studies, other authors state that freezing does not have any influence. For instance, Kiefer et al. [57], who performed indentation tests in bovine articular cartilage, concluded that the cryopreservation protocol used does not have any effect on its mechanical properties. Hongo et al. [58] obtained the same conclusion in porcine intervertebral discs, after conducting several types of test (lateral flexion, flexion-extension and rotation cycles) to analyse the effect of multiple freezing-thaw cycles. Van Haaren et al. [59] froze goat cortical bone during five different storage periods to study the effect on torsion, bending stiffness, strength and hardness. No significant differences were found. Nazarian et al. [60] carried out experimental tests on murine bone and stated that freezing does not change its mechanical properties. Szarko et al. [61] used indentation on bovine articular cartilage to check the influence of different storage temperatures. For −20◦C and −80◦C, they did not found alteration in the mechanical properties. Wex et al. [62] studied the mechanical properties of porcine liver tissue at different temperatures with shear relaxation and rheometry tests. No influence was noticed. Torimitsu et al. [63] froze cadaver skulls during different periods of time, concluding that it has little effect on the mechanical properties after conducting destructive 3-point bending tests. Wieding et al. [64], who performed destructive 4-point bending tests in ovine cortical bone using three different preservation techniques, did not find any difference in the elastic properties. Obviously, the best option is testing the specimens submerged in saline solution at 37◦C and as fresh as possible, trying to avoid freezing or any other preservation technique. In general, it is a good practice to avoid any procedure whose influence is not clear and could amplify the scatter of the results, but unfortunately, this is not always possible. 4.1.1.2 Zeroing the displacement Another important point to consider is the definition of the zero displacement in the test. A wrong definition of this zero can produce errors in the characterisation of the mechanical behaviour of the tissue. For example, in a uniaxial compression test between two platens, if the upper platen is the movable one, while bringing the platens together before the experiment it is crucial to know where it should be stopped, which defines the zero of the displacement. Errors in this estimation may lead to a prestress if the platen has been moved too much towards the specimen or a shift in the measured displacement if the platen has been stopped too far away from the specimen. This error can have a big influence on the results, mainly in the toe region. Usually, the problem lies in the fact that the top surface of the specimen is not flat, which complicates the visual detection of contact. Moreover, since the contact cannot be considered as fully established when the platen just touches the upper part of this top surface, the platen must be positioned such that it is compressing a certain portion of the specimen. This occurs, for example, in the articular discs of the temporomandibular joint, which has a sloping top surface. With this type of specimens, some authors measure their height in several points to calculate a mean height. Then, they zero the displacement when the distance between the 56 Mechanical properties of breast tissues platens is equal to the mean height [65]. In the specimens in which the top surface can be considered flat, the procedure is easier. Some authors move the upper platen until it comes into contact with the specimen, visually estimating this contact. Other authors lower the upper platen until it contacts with the specimen, determined from both visual inspection and from an initial increase in force (0.01N) [66]. In this thesis, the method to extract and prepare the specimens allowed obtaining them with the top and bottom surfaces plane and approximately parallel and, thus, for simplicity, the platen was positioned at the mean height of the specimen by visual inspection. 4.1.1.3 Scatter of the properties A significant point to take into account too is the scatter in the properties of biological tissues. The same type of tissue in different individuals may not have the same behaviour. Even the comparison of the tissue in different regions of the same individual can lead to different properties. The general structure of a given healthy tissue is the same in every person, but the specific composition (relative quantities of the inner components) can vary, depending on many factors: life style, diet, weight, etc. For example, the bone mineral density is highly dependent on the physical activity [67]. In the same person, the structure of a tissue can be quite different across regions. For instance, it is known that the structure and composition of the articular cartilage depends on the depth [68, 69]. This variation in properties between individuals and across regions introduces an additional difficulty to characterise a given tissue. This makes necessary the use of statistical tools to overcome this trouble, assigning average properties which are able to represent the whole tissue. 4.1.2 Computational estimation of properties Another way to characterise the behaviour of a biological material is a purely computational approach. It consists in knowing two different configurations of the body during the deformation process and simulating the transition from one to the other. To do that, it is necessary to choose a strain energy function beforehand and design an iterative process in which the material parameters of the strain energy function are changed in each iteration towards the optimisation. Finally, an optimum is reached if, starting from one configuration, the other configuration is reached (within a certain tolerance). In general, the experimental approach is preferred because it provides more realistic properties of the tissue, though all the problems pointed out in the previous paragraphs and also all the approximations and errors that are usually committed while carrying out experimental tests must be kept in mind. The computational approach is easier from a practical point of view. It also tends to be more phenomenological, because it consists in matching two deformed shapes by taking into account a global behaviour. Moreover, the properties experimentally obtained are valid, in general, for any other situation, while those obtained Introduction 57 through numerical approximation are only valid, commonly, in this particular application. 4.1.3 Mechanical properties of breast tissues. State of the art The material model chosen for the tissues is a key factor as said above, for its strong influence on the accuracy of the computational models [70, 71]. The breast is entirely composed of soft tissues which usually have a strong 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 should be added [72]. Moreover, depending on the internal microstructure, they can be considered as isotropic or anisotropic, for example, by modelling the collagen fibre directions [38]. Some information can be found in the literature about modelling the mechanical properties of fat and fibroglandular tissues. Most of the works modelled both tissues as isotropic and quasi-incompresible materials. Some considered the fat as anisotropic due to the connective tissue network inside the adipose tissue [7, 14], but no reliable work was found which experimentally obtained these anisotropic properties. Sarvazyan et al. [73] measured the Young’s moduli of fat tissue, fibroglandular tissue and tumorous carcinoma. Krouskop et al. [74] obtained experimentally the elastic modulus of the breast tissue subjected to compression loading. Wellman [75] studied the non-linear behaviour of fat and glandular tissue and used exponential curves to describe the stress-strain relationship. Lorenzen et al. [76] determined the linear elastic properties of breast tissues by means of elastography, including breast tumours. Samani et al. [24] considered them as elastic and determined experimentally the Young’s moduli of samples of both tissues. Some time later, Samani and Plewes [41] 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. [77, 51] and by O’Hagan and Samani [78, 79] have provided the properties of breast tissue with tumoral inclusions. Other authors have tried to obtain the mechanical properties through numerical methods. Usually, they acquire Computed Tomography (CT) or MR images from the body in two different positions, for example, supine and standing up. Then, a FE model is constructed, and the constants of the strain energy function are determined by fitting the relative displacements between both positions. P´erez del Palomar et al. [18] fitted the constant of a neo-Hookean model for the tissue inside the breast with this technique (they did not differentiate between fat and glandula, considering them as a whole tissue). Rajagopal et al. [16, 17] also considered a unique neo-Hookean material and obtained the value of its constant starting from an unloaded configuration and reaching the prone position. Han et al. [7] obtained the properties of fat, glandula, muscle and tumour matching images in a mammographic compression of the breast and using a simple anisotropic model. Some time later, Han et al. [15] carried out a similar analysis, but using a neo-Hookean law and applying gravity loads in the prone position to reach the supine one instead of a compression. Gamage et al. [19] also simulated a change of position from prone to supine fitting a neo-Hookean model too. 64 Mechanical properties of breast tissues Fat volume ratio C10 C01 C11 C20 C02 (10−4MPa) RMSE (kPa) f= 0 (pure gland) 3.3 2.8 44.9 77.2 94.5 (A.Samani, 2004) f= 0.1µ3.2835 2.8271 42.5265 72.9896 89.2642 3·10−4 σ0.0074 0.0054 0.1206 0.0632 0.0521 f= 0.3µ3.2453 2.8826 37.8507 64.6884 79.0541 9·10−4 σ0.0089 0.0064 0.1857 0.0832 0.0911 f= 0.5µ3.2010 2.9331 33.3466 56.5675 69.2128 7·10−4 σ0.0166 0.0138 0.2114 0.1170 0.0911 f= 0.7µ3.1738 2.9600 28.7921 48.8806 59.9891 8·10−4 σ0.0179 0.0136 0.2476 0.1341 0.1146 f= 0.9µ3.1249 2.9899 24.5546 41.5000 51.2965 3·10−4 σ0.0059 0.0042 0.0967 0.0475 0.0428 f= 1 (pure fat) 3.1 3.0 22.5 38.0 47.2 (A.Samani, 2004) Table 4.2: Mean, µ, and standard deviation, σ, obtained for each constant of the polynomial strain energy function. The fat volume ratio is termed as f. Figure 4.2: Correlations between C10 and C01 and fat volume ratio, f. 4.2.1.2 Results The constants of the polynomial strain energy function were fitted for each of the 16 cases studied per fat volume ratio, by minimizing the SE given in equation (4.1). The mean and standard deviation for these constants were calculated for each fat volume ratio and are shown in table 4.2. The root-mean-squared errors (RMSE) for the 16 cases are also shown in table 4.2. Linear or quadratic regressions were selectively proposed between the 5 constants Cij and f(see figures 4.2 and 4.3). The small errors found in the least squares fitting of the constants with the three load cases applied (see table 4.2) 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 [41]. The mean values of the constants Cij vary much with the fat content, confirming that the fat volume ratio Model for the mixture of fat and glandular tissue 65 Figure 4.3: Correlations between C11,C20 and C02 and fat volume ratio, f. 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. 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%. 4.2.2 Study B 4.2.2.1 Methods After estimating the constants of the homogenized material for each fat volume ratio, other material distributions and geometries were simulated to test the behaviour of the homogenized material. In study B, the same procedure followed in study A was 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 was analysed in this study. Two models, with 1000 and 125 elements, respectively, were tested and compared to the original model of 8000 elements (see figure 4.4). 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 66 Mechanical properties of breast tissues Figure 4.4: Cubic models with 125, 1000 and 8000 elements. Number of elements C10 C01 C11 C20 C02 (10−4MPa) 8000 3.2010 2.9331 33.3466 56.5675 69.2128 1000 3.2056 2.9325 33.3109 56.4507 69.0912 125 3.2007 2.9421 33.5065 56.4038 68.9541 Table 4.3: Comparison between the uniform random distribution models with 8000, 1000 and 125 elements for f= 0.5. little lactiferous ducts, which join those TDLUs, have an approximate diameter of 40-60 µm [27, 29]. During breastfeeding, these ducts increase their size, but that can be considered a marginal state of the female breast. 4.2.2.2 Results 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 4.3. The constants of the strain energy function are very similar in all cases, confirming that the homogenization is also valid in a breast in which larger pieces of a single tissue exist. 4.2.3 Study C 4.2.3.1 Methods In study C, a simplified geometry of the tissue inside the breast, consisting in a spherical cap of 8 cm in radius and a height of 6 cm 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 4.5 and has 45.3% of fat. To construct this ramification, the point N1 (see figure 4.5) was established as the origin of a number of beams equally spaced over the spherical cap. The glandular tissue Model for the mixture of fat and glandular tissue 67 Figure 4.5: Ramified model. was assigned to the elements within the beams, whose width was arbitrarily set, such that gland and fat were approximately divided in equal proportions.1This 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. 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 [18, 21]. Both models were subjected to gravity loads simulating three body positions: supine (gravity acting in the negative xdirection of figure 4.5), 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 [102, 103, 104]. For the homogenized material, the density was calculated using the rule of mixtures. 4.2.3.2 Results The results of the spherical cap geometry are presented here. The displacement of the nipple (node N1, see figure 4.5) in both models (ramified and homogenized) for the three body position simulated are compared in table 4.4. The relative difference between those displacements is also shown in the table. The relative differences in displacement between both models (ramified and homogeneous), shown in 1The 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. 68 Mechanical properties of breast tissues Figure 4.6: Boundary conditions of the back plane in study D. table 4.4, 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 [19]. 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. 4.2.4 Study D 4.2.4.1 Methods 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 the right part of figure 4.6, the nodes that are highlighted in red color), 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 producing 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.05 N/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. Model for the mixture of fat and glandular tissue 69 Homogenized model Ramified model Relative difference Spherical cap Supine 1.83 1.85 0.86% Prone 2.02 2.04 0.83% Standing up 5.98 6.03 0.79% Spherical cap with more flexible boundary conditions Supine 17.665 17.591 0.42% Prone 13.491 13.463 0.2% Standing up 7.546 7.611 0.85% Table 4.4: Modulus of the displacements (mm) of the nipple in the spherical cap geometry and in the spherical cap geometry with more flexible boundary conditions in both models and in the three body positions simulated. Relative difference between those displacements. 4.2.4.2 Results The results of the spherical cap geometry with more flexible boundary conditions are presented here. The displacements of the nipple in both models for the three simulated body positions are compared in table 4.4. The relative difference between those displacements is also shown in the table. 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 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. 4.2.5 Some remarks on the homogenized model The usefulness of this homogenization study 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, as has been stated in section 2.6. For some purposes, the aforementioned fat proportion categories (0-25%, 25-50%, 50-75%, 75-100%) 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. The study has also some limitations. 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 [41], which have been used to develop the present model. It is also important to note that the strain energy function used by these authors is not policonvex. However, this fact does not necessarily mean that this function should be ruled out, as explained in [105, 106, 107]. 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 70 Mechanical properties of breast tissues would fail to describe the behaviour under very localized loads, like a needle insertion for example. 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. 4.3 Experimental characterisation of breast and abdominal tissues Due to the scarce information that can be found in the literature about the mechanical properties of the adipose tissue in the breast as well as in the abdomen, experimental tests were carried out. Besides, these tests were motivated by the fact that, using the properties of fat provided in the literature to build a FE model of the breast, its behaviour is very far from the real behaviour of this organ, as will be seen later in the chapter devoted to numerical simulations. Indeed, the deformation seen in the model using the constants given by Samani and Plewes [41] is very low, indicating that the behaviour with such constants is much stiffer than the real one. Therefore, the suitability of the properties provided in the literature needs to be checked with care. Even less information exists about the properties of glandula, because samples are very difficult to be obtained. Due to the fact that the properties of glandula and fat are similar [41], that many authors model the breast with only one material, merging together adipose and glandular tissue [18, 21, 99, 100, 17], and that the homogenization procedure explained in section 4.2 seems to be adequate to model the mixture of adipose and glandular tissues, a unique material is going to be considered inside the breast. 4.3.1 Test protocol 4.3.1.1 Preparation of specimens The adipose tissue samples were extracted from the breast and the abdomen of several patients subjected to a mastectomy surgery. The abdominal tissue used in the experiments was the portion of the abdomen which was not eventually used for the reconstruction (see figure 4.7). Therefore, its size depended on how the surgery was performed. In fact, it is not easy to acquire it because in many occasions there is not enough remaining tissue, or even the piece is so small that only a few specimens can be extracted from it, being useless to carry out subsequent statistical analyses. The breast tissue was usually obtained from the contralateral breast, when the doctors decided to apply a mastopexy to improve the shape and symmetry of both the healthy and reconstructed breasts (see figure 4.8). Normally, it is difficult to obtain large pieces of abdominal and breast fat from the same patient. The tissue was extracted from a total of four patients, with a mean and standard deviation age of 52.3±6.9 years. The age of the patients ranged from 42 to 57 years. Once the pieces were extracted in the surgery, they were introduced in a cool-box with ice to preserve Experimental characterisation of breast and abdominal tissues 71 Figure 4.7: One of the pieces of abdominal fat from which specimens were extracted. Figure 4.8: Breast fat sample. them during transportation from the hospital to the Mechanical Engineering Laboratory. This was done shortly after excision (as soon as possible) to minimise the time the tissue was at room temperature. In the lab, the skin was removed from the piece by cutting a slice of tissue of approximately 5 mm in thickness. In the case of abdominal fat, the rest of the piece was divided in two slices of approximately equal size (superficial, underneath the skin; and deep, on top of the abdominal muscle), and each slice was, in turn, divided in two parts (medial and lateral). See figure 4.9 for a clear description of these divisions. The slices had a thickness ranging from 5 to 10 mm. This corresponded to the height of the tested specimens and was limited to that range to avoid slenderer specimens that could lead to buckling problems. This height is equal or lower than that used by Miller-Young et al. [66]. The division into lateral and medial parts was not performed in the smallest pieces. In the case of breast fat, only the division into superficial and deep slices could be done for size limitations. These divisions were intended to carry out statistical analysis to check if the mechanical properties of the fat are different depending on the region. All these cuts were made with a meat slicer machine and a sharp knife. It was important to make them while the tissue was cool, though 72 Mechanical properties of breast tissues Figure 4.9: Diagram of the divisions in the sample. Figure 4.10: Cylindrical specimen, top and lateral views. not completely frozen. If not, the fat is really soft and really difficult to cut in a proper way. The reconstruction surgery with autologous tissue is quite long, so that the pieces could not be collected, transported to the lab, cut and tested in the same day of extraction. They were tested the days after, as soon as possible to reduce the time elapsed from excision to test. To avoid the degradation of the tissue, it was frozen. Each slice of tissue was wrapped in saline-soaked gauze, enveloped in a plastic film and introduced in hermetic vials to prevent dehydration and finally frozen at −20◦C. The tests were carried out in cylindrical specimens that were extracted from the slices with a hollow punch of 19mm of diameter. This needed to be done with the tissue frozen. Otherwise, the final shape of the specimens were irregular and far from cylindrical, because the tissue was largely deformed by the punch. Then the specimens needed to be extracted while the slice was still frozen. After that, the slice was returned to the freezer to prevent it from thawing. It is advisable to drive the punch into the tissue slowly, pushing with the hands. If this is done quickly, for example by hitting the punch with a hammer, the shape of the specimen is worse. A specimen with the final cylindrical shape can be seen in figure 4.10. Next, the specimen was submerged in saline solution at room temperature and allowed to thaw. Then, it was photographed with a rule next to it (to have a reference measurement) to measure its area through imaging techniques. It is important to make the photograph with the rule positioned at the height of top face of the specimen, in order to measure the true size of the specimen without perspective artifacts. It is also a good practice to Experimental characterisation of breast and abdominal tissues 73 Figure 4.11: Scheme of the test. Remarked: (1) load cell, (2) sample, (3) upper platen, (4) temperature controller, (5) acquisition system. take the photograph with the camera as far as possible from the specimen, to minimise likely optic errors. 4.3.1.2 Experimental setup With the specimens prepared as explained before, the relaxation tests were carried out. The specimens were compressed between two platens, in a uniaxial compression test. To do this, a servo-hydraulic testing machine (858 Mini Bionix II, MTS, Eden Prairie, USA) was used together with ancillary equipment specially designed to perform the test under certain conditions of humidity and temperature. A scheme of the experimental setup can be seen in figure 4.11. The experimental tests were carried out as follows. First, both platens of the testing machine were brought into contact to zero the displacement. Then, vaseline was spread on their surfaces to reduce friction (although this is probably not needed as the fat seemed to lubricate the platens itself). Once the specimen was thawed, it was placed in the center of the inferior platen of the testing machine. The specimen was not glued to the platen unlike in similar tests performed on articular discs in a previous study [65]. This was done to avoid any undesired effect of the boundary conditions. In all the conducted tests, the specimens remained within the platens, so that it was not necessary to glue them to prevent slipping. Next, a recipient was filled with saline solution (0.9% w/v of NaCl) at 37 ±1◦C and the specimen remained submerged in it during the whole test. The temperature was controlled with a heater and a thermostat. After having zeroed the force, the upper platen was manually lowered, slowly approaching the sample, and visually positioned in contact with the top surface of the sample, as made in [66]. That served to measure the thickness of the specimen and to find the starting point of the test. A preconditioning was applied to each sample: 20 cycles from 0% to 10% strain at 1 Hz, like in [108]. The preconditioning was followed by a ramp from 0% to strain, where was 40%, 50% or 60% depending on 80 Mechanical properties of breast tissues + m X j=1 e−∆t/τjJ−1Fn+1Qj nFT n+1 +e−∆t/2τjβ∞ j(σ∞ iso n+1 −J−1Fn+1S∞ iso nFT n+1).(4.12) 4.3.4 Description of the experiments Following the test protocol shown in section 4.3.1 and using the fitting algorithm proposed in sections 4.3.2 and 4.3.3, the viscoelastic properties of the abdominal and breast fat were determined. The information about the patients, the samples and the test conditions is summarised in table 4.5. First, it is important to check which one of the four proposed strain energy functions fits best the experimental curves for each viscoelastic model. To do this, the mean CV was calculated for each group of specimens (sample division in table 4.5) and each strain energy function used. Then, the mean CV was calculated for each of the four strain energy functions and they were compared to determine which was the lowest CV and, therefore, the best model. The validity of both, the QLV and the IVV models, must be verified once the most suitable model for the elastic response is selected. This must be done by checking if the material parameters are independent of the strain rate and the strain level of the relaxation test. To carry out this validation, in a given piece of tissue from the same patient, three specimen groups were extracted and each one was tested with a different strain rate. Then, in a piece of another patient, three more groups were extracted and each one tested with a different strain level. Choosing specimens from the same patient avoids the problem of inter-individuals differences, but increases the difficulty to obtain a suitable piece of tissue, due to the fact that its size must be large enough to host the sufficient number of specimens needed for statistical purposes. The strain rates selected were 50% /s, 60% /s and 70% /s and the strain levels, 40%, 50% and 60%. If no difference exists between the material parameters obtained for the different strain levels and strain rates, that will mean that the model is valid to represent the behaviour of the adipose tissue. These two comparisons were carried out with statistical methods to compare the differences between groups. As there are multiple dependent variables but only one independent variable with several levels (the strain rate in one case and the strain level in the other comparison), the first choice was to perform a one-way multivariate analysis of variance (MANOVA) test. Of course, all the assumptions of MANOVA tests should be accomplished. If not, a non parametric MANOVA should be used instead. After checking the validity of QLV and IVV models, they were compared in terms of goodness of fit. The one able to better reproduce the viscoelastic behaviour of the adipose tissue was chosen as the final model. Finally, using this model, some more comparisons were carried out to check if the adipose tissue has the same properties in different regions of the abdomen, and to investigate the differences between the adipose tissue of the abdomen and the breast. Experimental characterisation of breast and abdominal tissues 81 Patient Age Extraction area Sample division Test strain rate Test strain level Noof specimens A 42 Abdomen Superficial 50% /s 50% 8 Deep 50% /s 50% 10 B 55 Abdomen Superficial-lateral 50% /s 50% 14 Superficial-medial 50% /s 50% 26 Deep-lateral 50% /s 50% 17 Deep-medial 50% /s 50% 23 Breast Superficial 50% /s 50% 12 Deep 50% /s 50% 14 C 55 Abdomen Superficial-lateral 60% /s 40% 11 50% 11 60% 11 Superficial-medial 60% /s 40% 12 50% 12 60% 10 Deep-lateral 60% /s 40% 11 50% 11 60% 10 Deep-medial 60% /s 40% 10 50% 10 60% 9 D 57 Abdomen Superficial-lateral 50% /s 50% 13 60% /s 12 70% /s 12 Superficial-medial 50% /s 50% 12 60% /s 12 70% /s 12 Deep-lateral 50% /s 50% 11 60% /s 11 70% /s 11 Deep-medial 50% /s 50% 11 60% /s 10 70% /s 11 Table 4.5: Information of patients and tissue extracted. 82 Mechanical properties of breast tissues 012345678 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (s) Cauchy stress (kPa) Strain rate: 60% /s. Strain level: 40% Experimental stress (˜σ) Analytical stress (isotropic GOH) Analytical stress (Exponential) Analytical stress (Ogden) Analytical stress (5 terms polynomial) Figure 4.15: Example of an experimental stress record fitted with different hyperelastic models, for the QLV model. 4.3.5 Results The results obtained in the experimental tests of abdominal and breast fat are presented in this section. Each specimen was tested and the results treated with the algorithm, obtaining the material parameters of the hyperelastic model used for Teand the constants giof the Prony series in the case of the QLV model; and the material parameters of the hyperelastic model and the constants β∞ jin the case of the viscoelastic model with internal variables. 4.3.5.1 Goodness of fit of the hyperelastic models with the QLV model Figure 4.15 compares a typical experimental stress record, ˜σ, with the fitting curves for each one of the strain energy functions proposed. In figure 4.16 a detail of the loading ramp can be seen. In these figures, it can be observed that the toe region of the experimental curve was not very wide in comparison with other tissues, as for example the articular discs of the temporomandibular joint [65]. The loading ramp applied to the tissue was quickly noticed by an increase in its stress state. The relaxation was very quick, with a high relaxation percentage occurring few seconds after the peak. Although it is not shown here, it could be observed that more than a 90% of the peak stress was relaxed after 15 minutes. Comparing these stress relaxation curves with those obtained by Miller-Young et al. [66] testing specimens of the human calcaneal fat pad in compression, the peak stress of the present work was of the same order of magnitude and in both, 75% of the stress was relaxed within the first minute. The fitting was very accurate for all cases and quite similar for all the fitted hyperelastic models. In both figures 4.15 and 4.16, it can be noticed that the isotropic GOH and the 5 terms polynomial models Experimental characterisation of breast and abdominal tissues 83 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (s) Cauchy stress (kPa) Strain rate: 60% /s. Strain level: 40% Experimental stress (˜σ) Analytical stress (isotropic GOH) Analytical stress (Exponential) Analytical stress (Ogden) Analytical stress (5 terms polynomial) Figure 4.16: Detail of figure 4.15 to notice the loading ramp. Patient Extraction area Model Whole curve CV (%) Loading ramp CV (%) A Abdomen Exponential 4.13 18.27 Ogden 4.15 18.09 GOH 3.21 11.82 Polynomial 3.05 10.89 B Abdomen Exponential 4.48 17.22 Ogden 4.31 15.90 GOH 3.26 8.58 Polynomial 3.15 7.91 Breast Exponential 3.44 12.85 Ogden 3.35 12.06 GOH 2.65 6.83 Polynomial 2.56 6.23 C Abdomen Exponential 4.95 21.46 Ogden 4.96 21.22 GOH 4.11 15.38 Polynomial 3.86 14.13 D Abdomen Exponential 4.16 19.88 Ogden 4.05 18.70 GOH 3.06 11.49 Polynomial 2.92 10.67 Table 4.6: Coefficient of variation for each patient and model, for the QLV model. 84 Mechanical properties of breast tissues are closer to the experimental curve than the Ogden and exponential ones, although not too much. The average coefficient of variation, CV , for each patient and model is presented in table 4.6. The CV has been evaluated in the whole curve and in the loading ramp interval separately. In view of table 4.6, it is clear that the best model is the polynomial one, closely followed by the GOH model, and a little bit farther by the Ogden and exponential models, in this order. However, the CV is very low for each patient and model. This fact means that the fitting was quite good in all cases. However, due to the fact that the polynomial model has so many constants, in the results given by the fitting algorithm several constants from the total of five were zero for the majority of the specimens. Moreover, the constants equal to zero were not the same for all of the specimens. This means that a different model is being used to fit each specimen. For example, if constants C10,C11 and C02 are zero, the strain energy function which is being fitted is: Ψ = C01(I2−3) + C20(I2−3)2,(4.13) whereas if constants C01 and C20 are zero, the strain energy function which is being fitted is: Ψ = C10(I1−3) + C11(I1−3)(I2−3) + C02(I1−3)2.(4.14) Depending on the particular specimen, some constants are zero but not the same ones in all specimens. Then, if an average is calculated for each constant, a non zero mean value results. As a consequence, the whole material is represented by a 5 terms polynomial. But if the curve calculated with the mean constants is compared with the curves of each individual specimen, big differences can be found. This is because all the specimens have been fitted with a different strain energy function. From a physical point of view that does not make much sense. A specific material should be represented by a unique function. Also, neither equation (4.13) nor (4.14) correspond to well-known or widely used strain energy functions, which does not help to the choice of the 5 terms polynomial function as a suitable model for the fat. With the GOH model, something similar occurs. In many specimens, the fitted constant k2was equal to zero. The main problem of k2= 0 is that it generates an indeterminate form (0/0) in the strain energy function. If the limit of the function when k2→0 is taken: lim k2→0 k1 2k2hek2(I1−3)2−1i=k1 2(I1−3)2.(4.15) Therefore, the GOH strain energy function (equation (4.2c)) would be in this case: Ψ = C10(I1−3) + k1 2(I1−3)2,(4.16) Experimental characterisation of breast and abdominal tissues 85 Strain level µ(kPa) α g1g2g3g4g5g∞ 40% 2.810 11.919 0.509 0.342 0.065 0.028 0.028 0.028 ±2.608 ±1.457 ±0.262 ±0.193 ±0.033 ±0.015 ±0.016 ±0.017 50% 8.757 7.280 0.606 0.290 0.045 0.019 0.021 0.018 ±4.874 ±1.229 ±0.198 ±0.175 ±0.017 ±0.007 ±0.009 ±0.008 60% 18.860 5.048 0.659 0.249 0.040 0.018 0.019 0.014 ±16.274 ±0.964 ±0.175 ±0.145 ±0.016 ±0.007 ±0.007 ±0.006 Table 4.7: Mean ±standard deviation of the QLV constants for the different groups, for patient C. Strain rate µ(kPa) α g1g2g3g4g5g∞ 50% /s 6.277 8.504 0.579 0.290 0.057 0.027 0.025 0.020 ±4.076 ±1.750 ±0.215 ±0.170 ±0.027 ±0.014 ±0.013 ±0.011 60% /s 8.536 8.316 0.610 0.264 0.054 0.026 0.025 0.020 ±14.546 ±1.953 ±0.193 ±0.139 ±0.027 ±0.014 ±0.013 ±0.011 70% /s 6.317 8.564 0.606 0.263 0.055 0.029 0.027 0.020 ±4.343 ±1.683 ±0.173 ±0.141 ±0.022 ±0.014 ±0.011 ±0.008 Table 4.8: Mean ±standard deviation of the QLV constants for the different groups, for patient D. which is a particular case of the 5 terms polynomial strain energy function. This is the reason why the GOH and polynomial models perform very similarly when fitting the experimental curves. Therefore, the isotropic GOH model is not appropriate to model the fat either. In conclusion, as some problems were found with these two models, the Ogden model was selected to represent the behaviour of the adipose tissue. 4.3.5.2 Validity of the QLV model Once a suitable model has been selected to represent the behaviour of the adipose tissue, it has been used to check the validity of the QLV approach. As stated before, this must be done by checking if the fitted model constants are independent of the strain level and the strain rate. For that, the specimens extracted from patient C (to check the independence of the strain level) and those extracted from patient D (to check the independence of the strain rate) were used. The validity was proved only in the abdominal tissue, because it was difficult to obtain a large piece of breast fat to have enough specimens to test. Extracting the specimens from several patients was considered inadequate because differences between individuals could exist. In table 4.7, for patient C, the mean and standard deviation for each strain level and each constant are presented. For patient D, the same is done in table 4.8. In both tables, it can be noticed the high values of the constants g1and g2, which correspond to the relaxation times τ1= 0.01 and τ2= 0.1. These high values highlight the fast relaxation in the short-term. To check if the material constants are independent of the strain level and the strain rate, a MANOVA 86 Mechanical properties of breast tissues test was performed for each patient. In the specimens extracted from patient D, the independent variable was the strain rate with three levels: 50% /s, 60% /s and 70% /s; and the dependent continuous variables (DVs) were the seven QLV constants: µ,α,g1,g2,g3,g4and g5. The material parameter g∞was not included in the statistical analysis because it is a linear combination of the rest of constants gi, due to the normalization condition of G(t): G(0) = 1 = g∞+ 5 X i=1 gi.(4.17) Before applying a MANOVA, some assumptions needed to be checked. First, multinormality was checked for each independent group using the test developed by Cardoso de Oliveira and Ferreira [116]. This test was not significant for the 50% /s group (p=.098) and 60% /s group (p=.160), but significant for the 70% /s group (p < .001) and, thus, multinormality was violated, making it necessary to perform a non-parametric MANOVA (NMANOVA). The seven variables were compared in the NMANOVA: µ,α,g1,g2,g3,g4and g5. The NMANOVA test performed in this thesis is a multivariate extension of the Kruskal-Wallis test, developed by Katz and McSweeney [117]. No significant differences were found for the three groups compared (p=.819). Some authors state that MANOVA is indicated if the dependent variables are correlated, but not so strongly that multicollinearity exists [118, 119]. As this is not the case for all the variables in this study, the statistical analysis was repeated after eliminating the correlated variables, just to be sure of the correctness of the previous conclusion. Constants g1and g2were strongly correlated and the same occurred to g3,g4and g5(Spearman |R|> .86 in both cases). Thus, they could be regarded as the same variable for statistical purposes. Then, the analysis was repeated with 4 DVs: µ,α,g1and g3. The conclusion was the same: there were no statistical differences. Therefore, the material constants could be considered independent of the strain rate. In the specimens extracted from patient C, the independent variable was the strain level with three levels: 40%, 50% and 60%; and the dependent continuous variables (DVs) were the seven QLV constants: µ,α,g1, g2,g3,g4and g5. The material parameter g∞was not included either. Multinormality was checked for each independent group. This test was significant for the three groups: 40% (p < .001), 50% (p=.018) and 60% (p=.020) and, thus, multinormality was violated, making it necessary to perform a NMANOVA. The seven variables were compared in the NMANOVA: µ,α,g1,g2,g3,g4and g5with the test developed by Katz and McSweeney [117]. Significant differences were found (p<.001). Afterwards, post-hoc tests were carried out to detect the origin of these differences. To do so, the proposal by Katz and McSweeney [117] was followed again. Significant differences were found for µ(between 40% and the other two groups, p<.001), α(between all groups, p<.001), g1(between 40% and 60%, p=.038), g3(between 40% and Experimental characterisation of breast and abdominal tissues 87 012345678 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (s) Cauchy stress (kPa) Strain rate: 60% /s. Strain level: 40% Experimental stress (˜σ) Analytical stress (isotropic GOH) Analytical stress (Exponential) Analytical stress (Ogden) Analytical stress (5 terms polynomial) Figure 4.17: Example of an experimental stress record fitted with different hyperelastic models, for the IVV model. 60%, p < .001) and g4(between 40% and 60%, p=.018). As in the previous case, the statistical analysis was repeated after eliminating the correlated variables, just to be sure of the correctness of the previous conclusion. In this case, µand αwere strongly correlated, and the same occurred to g1and g2; and to g3, g4and g5(Spearman |R|> .80 in all cases). Therefore, the analysis was repeated with 3 DVs: µ,g1and g3. The conclusion was the same: there were statistical differences (p < .001). As a consequence, the material constants are not independent of the strain level. With these statistical results, it can be concluded that the QLV model is not suitable to model the adipose tissue, due to the fact that the independence of the model with the strain level cannot be ensured. Moreover, the assumption of the multiplicative decomposition of the stress into a function which only depends on time and another one which only depends on strain (equation (3.59)) is not valid, due to the fact that ¯ G(t), by definition, can not depend on the strain level. 4.3.5.3 Goodness of fit of the hyperelastic models with the IVV model Figure 4.17 compares a typical experimental stress record, ˜σ, with the fitting curves for each one of the strain energy functions proposed. In figure 4.18 a detail of the loading ramp can be seen. The experimental curve shown in the figures is the same one shown in figures 4.15 and 4.16. As in the QLV model, it can be seen that the fitting is very accurate for all cases and quite similar for all the fitted hyperelastic models. In both figures 4.17 and 4.18, it can be noticed that the isotropic GOH and the 5 terms polynomial models are closer to the experimental curve than the Ogden and exponential ones, although not too much. The goodness of fit is slightly better than with the QLV model for all the 88 Mechanical properties of breast tissues 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 −8 −7 −6 −5 −4 −3 −2 −1 0 Time (s) Cauchy stress (kPa) Strain rate: 60% /s. Strain level: 40% Experimental stress (˜σ) Analytical stress (isotropic GOH) Analytical stress (Exponential) Analytical stress (Ogden) Analytical stress (5 terms polynomial) Figure 4.18: Detail of figure 4.17 to appreciate the loading ramp. hyperelastic strain energy functions tried, except for the polynomial model, in which case it is slightly worse, although that is difficult to notice with the naked eye. It can be noticed by looking at the average coefficient of variation, CV , which is presented for each patient and model in table 4.6, for the QLV model, and in table 4.9, for the IVV model. The CV has been evaluated again in the whole curve and in the loading ramp interval separately. Looking at table 4.9, it is clear that the best model is the GOH model, closely followed by the polynomial one, and a little bit farther, by the Ogden and exponential models, in this order. However, the CV is very low for each patient and model. This fact means that the fitting was quite good in all cases. As in the previous QLV model, the use of the polynomial and the GOH strain energy functions lead to many constants equal to zero, causing the same problems discussed in section 4.3.5.1. Therefore, the Ogden model was selected to represent the behaviour of the adipose tissue, as done with the QLV model. 4.3.5.4 Validity of the IVV model Once a suitable model has been selected to represent the behaviour of the adipose tissue, it has been used to check the validity of the IVV approach. As stated before, this must be done by checking if the fitted model constants are independent of the strain level and the strain rate. For that, the specimens extracted from patient C and those extracted from patient D were used, as in the validation of the QLV model. In table 4.10, for patient C, the mean and standard deviation for each strain level and each constant are presented. For patient D, the same is done in table 4.11. It can be noticed how the values of the variables β∞ jdecrease as jincreases. That is to say, the highest β∞ jcorrespond to the lowest relaxation times τj. Experimental characterisation of breast and abdominal tissues 89 Patient Extraction area Model Whole curve CV (%) Loading ramp CV (%) A Abdomen Exponential 3.79 16.30 Ogden 3.86 16.39 GOH 3.16 11.54 Polynomial 3.18 11.88 B Abdomen Exponential 4.12 15.00 Ogden 4.05 14.23 GOH 3.27 8.65 Polynomial 3.38 9.58 Breast Exponential 3.18 10.73 Ogden 3.17 10.47 GOH 2.70 6.78 Polynomial 2.72 6.95 C Abdomen Exponential 4.50 18.84 Ogden 4.56 18.95 GOH 3.98 14.69 Polynomial 4.01 14.93 D Abdomen Exponential 3.81 17.66 Ogden 3.77 16.99 GOH 3.04 11.36 Polynomial 3.14 12.27 Table 4.9: Coefficient of variation for each patient, for the IVV model. Strain level µ(kPa) α β∞ 1β∞ 2β∞ 3β∞ 4β∞ 5 40% 0.063 11.481 93.284 12.269 3.632 1.601 1.597 ±0.040 ±1.676 ±57.32 ±5.702 ±1.078 ±0.502 ±0.472 50% 0.190 6.632 103.973 13.002 3.579 1.635 1.800 ±0.125 ±1.369 ±73.506 ±6.436 ±1.161 ±0.569 ±0.610 60% 0.545 4.381 96.235 13.242 3.820 1.938 2.030 ±1.316 ±1.421 ±51.204 ±6.136 ±0.776 ±0.454 ±0.349 Table 4.10: Mean ±standard deviation of the IVV constants for the different groups, for patient C. Strain rate µ(kPa) α β∞ 1β∞ 2β∞ 3β∞ 4β∞ 5 50% /s 0.121 8.165 84.956 15.444 4.278 2.151 1.980 ±0.068 ±2.006 ±61.953 ±7.511 ±1.479 ±0.715 ±0.455 60% /s 0.156 8.111 82.744 15.055 4.096 2.077 2.008 ±0.198 ±2.160 ±59.015 ±7.816 ±1.195 ±0.556 ±0.410 70% /s 0.130 8.249 83.547 13.209 4.159 2.281 2.114 ±0.069 ±1.893 ±81.419 ±6.727 ±1.602 ±1.058 ±0.544 Table 4.11: Mean ±standard deviation of the IVV constants for the different groups, for patient D. 96 Mechanical properties of breast tissues Superficial breast Deep breast Superficiallateral abdomen Superficialmedial abdomen Deep-lateral abdomen Deep-medial abdomen Superficial breast p=.111 p= 1 p=.001 p=.004 p < .001 Deep breast p= 1 p= 1 p= 1 p= 1 Superficiallateral abdomen p=.069 p=.032 p=.015 Superficialmedial abdomen p=.032 p=.192 Deep-lateral abdomen p= 1 Deep-medial abdomen Table 4.14: P-values of the post-hoc comparisons, showing the lowest p-value of the seven dependent variables in each cell and highlighting in bold typeface the significant ones (p<.05). Also important although with less clinical relevance are the differences found between regions of the abdominal adipose tissue. It seems that no differences exist in the mechanical properties between both regions of the superficial layer, between both parts of the deep layer and between both parts of the medial regions. However, differences between the deep and superficial layers were founded, that is to say, the mechanical properties of the abdominal adipose tissue seem to change with the depth. This fact is in accordance with other authors’ findings [47] who suspected it, although they did not proved it. In addition, all the differences were detected in the constants which control the stress relaxation. Apparently the elastic constants were the same for the whole adipose tissue in the abdomen and breast, but the viscous constants were not. This implies that the behaviour under static loading can be consider equivalent, but not under dynamic loading. In figure 4.20, as previously, the lowest and the highest experimental curves, in terms of the peak Cauchy stress, are drawn for the superficial breast fat and the three abdominal groups with which differences were found. The difference between the superficial breast fat and both deep regions of the abdomen is clear in view of the different ranges they have. The difference between the superficial breast and the superficial-medial abdomen is not easily seen with the naked eye. For practical reasons, the whole set of experimental curves could not be shown, but their general trend is as follows. Most of the specimens of the superficial-medial abdomen lie in the region enclosed by the curves of the deep layers of the abdomen. Moreover, most of the specimens of the superficial breast fat lie closer to its lowest limit than to the highest one, which could explain the origin of the differences found with the statistical analysis. In figure 4.21, a figure similar to the previous one is presented, but now for the four abdominal groups. The differences between both deep groups and the superficial-lateral one seems to be justified by the differences Experimental characterisation of breast and abdominal tissues 97 0 10 20 30 40 50 −30 −25 −20 −15 −10 −5 0 Time (s) Cauchy stress (kPa) Superficial breast Superficial-medial abdomen Deep-lateral abdomen Deep-medial abdomen Figure 4.20: Comparison between the experimental curves for the superficial breast and the three abdominal groups with which differences were found. 0 10 20 30 40 50 −30 −25 −20 −15 −10 −5 0 Time (s) Cauchy stress (kPa) Superficial-lateral abdomen Superficial-medial abdomen Deep-lateral abdomen Deep-medial abdomen Figure 4.21: Comparison between the experimental curves for the four abdominal groups. 98 Mechanical properties of breast tissues in the ranges. For the difference between the deep-lateral and superficial medial regions, although the ranges seem different, it is important to recall from the analysis of the previous figure that most of the specimens for the superficial-medial part lie in the region enclosed by the curves of the deep layers of abdominal fat. Therefore, no difference seems to be noticeable with the naked eye. If the curves for the specimens of both groups are studied in detail, it seems that they do not follow the same tendency in some areas. This could mean that the relaxation was different during the experimental test. It is important to note that this study should be broaden to a larger number of patients, to increase the certainty of the above conclusions. 4.3.6 Comparison between the Samani’s constants and the ones obtained here In this section, the constants given by Samani for the breast adipose tissue [41], which are the unique ones for breast fat in the literature, are compared with those obtained in this work. Here, the experimental tests were carried out on human adipose tissue mainly because of two reasons. First, because little information exists in the literature about the properties of this tissue, and even less about the breast and abdominal fat in particular. But also due to the fact that the FE models that used the experimental constants provided by Samani resulted very stiff compared with reality, as will be seen later in this document (chapter 5). When gravity loads were applied, the model was almost not deformed. That was actually the reason to carry out the previous experimental tests: to be sure of the goodness of the constants provided by Samani. Here, the viscoelastic properties of the adipose tissue have been obtained. To be able to do a comparison with the information given by Samani, only the elastic part has been taken into account. Of course, the constants compared in this section are those corresponding to the fitting of the IVV model with the Ogden strain energy function for the elastic behaviour, as this was the finally selected model. Moreover, as Samani did experiments in the breast tissue, only our results for the breast have been compared. In figures 4.22 and 4.23, the stress-stretch curve of a uniaxial test simulated in a material with a behaviour described by the polynomial model and Samani’s constants is compared with that simulated in a material having the behaviour fitted in the previous sections, both for the superficial and deep layers of the breast respectively. As can be seen in both figures, the difference between both models is very noticeable, even for small stretches. The behaviour with the Samani’s model is much stiffer. As will be seen later, after noticing that the FE model behaved much stiffer than in the reality, another set of simulations were conducted by dividing Samani’s constants by some factors (10,100...) to analyse the effect of the softening of the tissues. Then, another curve is added to figures 4.22 and 4.23 in which the same curve is plotted but with Samani’s constants divided by 100. As can be seen, the latter curve is very similar to that using the constants obtained in this thesis. In section 5.6.6, the constants will be compared Experimental characterisation of breast and abdominal tissues 99 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 −5 −4 −3 −2 −1 0 1 2 3 4 5 Stretch (λ) Cauchy stress (kPa) Samani’s constant Samani’s constant divided by 100 Constants of this work Figure 4.22: Comparison between the elastic constants given by Samani for the breast fat and the ones calculated in this work for the superficial layer. 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 −5 −4 −3 −2 −1 0 1 2 3 4 5 Stretch (λ) Cauchy stress (kPa) Samani’s constant Samani’s constant divided by 100 Constants of this work Figure 4.23: Comparison between the elastic constants given by Samani for the breast fat and the ones calculated in this work for the deep layer. 100 Mechanical properties of breast tissues in a FE model of the breast. Therefore, and based on our experimental results and on the computational work, it seems that the constants provided by Samani do not match very well to the real behaviour of the breast adipose tissue. Chapter 5 Finite element models of the breast 5.1 Introduction The cooperation between medicine and biomechanics has emerged as a major issue which can be very useful to widen the knowledge and drive the research. In particular, the FE models of biological systems have been widely used for many years and have proved to be very useful in many applications. Although it is quite complicated to model the behaviour of biological tissues, their interactions and functions, many advances have been carried out and a lot of help has been provided to clinicians in many fields. As stated in the introduction of this thesis, the motivation of this work is to shed light on the modelling of the behaviour of the healthy breast, because, to the best of the author’s knowledge, nobody has been able to simulate the natural movement of the breast, i.e., its global behaviour. Many breast FE models exist for different applications. Mainly, they are devoted to simulate mammographic compression, gravity loading or to match different medical images, as mammography and MRI. Azar et al. [3] constructed a FE model to predict the location of a tumour when the breast is compressed between two plates. They tested the model comparing the estimated location of the tumour with real MR images. In the same year, Samani et al. [5] made a model, also for compression between two plates, validating it with a phantom and finally with a real breast with the aid of MRI. Schnabel et al. [8] presented a method to validate nonrigid image registration using FE models. For many years, Ruiter el al. [9, 10] and Hopp et al. [11, 12, 13] have been working to develop medical image registration algorithms, using FE models of the breast to simulate breast mammographic compression. Pathmanathan et al. [6] built a breast FE model and obtained an undeformed configuration and a breast compression, but they did not compare the numerical results with experimental ones. Rajagopal et al. [16, 17] developed a subject-specific breast model subjected to gravity, using cubic-Hermite basis functions to generate the geometry. They identified the reference state as the breast immersed in water. Considering the breast composed of only one material, and modelling it 102 Finite element models of the breast with a neo-Hookean strain energy function, they fitted its material parameter to simulate the prone position. P´erez del Palomar et al. [18] constructed a FE model to simulate the change in position between supine and standing up. They measured the distances between some characteristic points in the model and compared them with the real ones. Lapuebla-Ferri [21] presented a model to simulate a prosthesis insertion for augmentation mammoplasty. The deformed shape of the FE model was compared with the real shape. Han et al [7] built a FE model for mammary compression between two plates and firstly compared different solvers. Finally, they compared the results between four different material characteristics: isotropic-anisotropic and homogeneous-heterogeneous. Gamage et al. [19] built a model to simulate the position change from prone to supine with heterogeneous FE models. They compared the results with the homogeneous models. Zain-Ul- Abdein et al. [20] constructed a FE model, starting from a reference position, that was the breast immersed in water. Then, they applied gravity loads to simulate the upright position and compared it with a real image through the distance between some points. Mertzanidou et al. [14] presented a framework for MRI to X-ray mammography registration, i.e., using a compression simulation. They optimised both material and pose parameters simultaneously during registration. Han et al. [15] proposed a hybrid nonlinear biomechanical model based registration method. It combined a biomechanical FE modelling with image registration techniques. The FE model gave a plausible deformation which was used to initialize the registration. Pianigiani et al. [4] developed a subject-specific breast model based on anthropometric measurements. They validated the model by comparing it with reality in a craniocaudal compression. Then, they simulated 3 different loading conditions, used 3 different materials and compared the displacement of a tumour inclusion between them. Some authors have worked trying to determine the reference state of the breast, i.e., the geometry of the breast if it were not subjected to any load. Rajagopal et al. [22] proposed an inverse method to obtain this unloaded configuration from a deformed state. They validated this method comparing the computed reference state in two problems: the analytic solution of pressure inflation and axial extension of a thickwalled cylinder and the experimental results of a silicon gel sample subjected to gravity loads. In the same year, Rajagopal et al. [16] identified the reference state as that where the breast is immersed in water. However, the density of the breast tissues is not equal to the water density and therefore the gravity is not completely cancelled. Moreover, their assumption does not take into account the possible prestress of the skin and breast tissues. Zain-Ul-Abdein et al. [20] considered also this situation as the reference state. Vavourakis et al. [23] recently proposed a displacement-pressure formulation to obtain the unloaded configuration of soft tissues and applied it to the breast. Other authors, such as Riveros et al. [120], have published methods to obtain the reference state which can be applied to the breast, though they were initially applied to other tissues. All these models can be used for the application they were designed for (with better or worse results FE model construction 103 depending on the particular model), but when they are used to simulate a different scenario, they normally do not work. P´erez del Palomar [18] and Zain-Ul-Abdein et al. [20] validated their models measuring the distance between some points of the mesh. But simply looking at the deformed shape of the models with the naked eye, it is clear that they are far from the real deformation. Many authors validated their models using a compression test and comparing the deformed shapes with a real image [3, 24, 7, 14, 4], but these deformations were obtained by imposing displacements in the contact between the breast and the plates and considering the materials as incompressible. Therefore, it was quite foreseeable that the deformed configurations were similar to reality, but not necessarily the stresses. As stated previously, our aim was to simulate the change of position from supine to prone and estimate the deformed shape of the breast due to the action of gravity loads. The initial objective was to simulate the behaviour under any change of position of the body, but to simplify the problem, only the supine-prone change was chosen because it is considered as the most complex, due to the great deformations which are produced. For that purpose, a FE model of a female breast was constructed and its deformed shape after applying gravity was compared with superficial 3D scan images taken from the same individual that the FE model. Using the material properties and boundary conditions provided in the literature, the behaviour of the breast was quite far from the real one. That could mean that those properties and boundary conditions existing in the literature might be useful for a limited number of applications, but not to reproduce the real breast behaviour in a large deformations problem like this. Therefore, this thesis intends to provide more information to comprehend the way the breast deforms under loading. Once the FE model was built, simple boundary conditions and material models were initially tried. Then, they were changed to study their influence and to approach a more realistic deformation. This chapter is structured so that the different models are presented following the real timeline, explaining the improvement provided by each model, and why the next model was built. First, it will be explained how the FE model was constructed and the possible boundary conditions applied as well as some other information needed to understand the whole process. 5.2 FE model construction The FE model was constructed from CT images taken from a 48 years old female patient at the Virgen del Roc´ıo Hospital before surgery for breast reconstruction. The resolution of these images is worse than that of MR images, so hindering the differentiation of soft tissues. However, the protocol of this hospital dictates the use of that imaging technique, which was therefore the only one available. The CT images were taken in supine position, so that this was the initial position of the FE model. The CT study consisted of 120 images with a resolution of 0.87 x 0.87 x 5 mm. The patient stretched her arms over her head (aligned with the trunk) and kept them extended during acquisition. A sample image of 104 Finite element models of the breast Figure 5.1: Sample image of the CT with the original resolution (0.87 x 0.87 mm in the transverse plane) viewed in craniocaudal direction from bottom to top. the original CT can be seen in figure 5.1. In it, the spine can be seen at the bottom of the image, the heart in the center-right part and the right breast is observed on the left. The left breast is not present because the CT was taken after a mastectomy surgery of this breast. The Simpleware software was used to process the images in order to build the FE model. In particular, the ScanIP and ScanFE modules were used. With the first one, the CT images were processed and the tissues segmented. Then, the file was exported to ScanFE, where the parameters to build the FE model were set, and finally the input file (.inp) generated for the simulation software Abaqus. The process is detailed next. First, the CT images were imported with the ScanIP module. This software works with voxels, whose dimensions correspond to the resolution of the images and from which the elements of the FE model were generated. As the initial resolution was 0.87 x 0.87 x 5 mm, the voxels had one dimension considerably larger than the others. To improve the aspect ratio of the voxels, the images needed to be resampled to convert the voxels into cubes. Another issue to take into account was that the smaller the voxels were, the higher the number of elements would be, so increasing the computation time. Therefore, it was important to reach a compromise between computation time and resolution of the model. A mesh convergence analysis was carried out to determine a mesh density which accomplished this compromise solution. At first, the images were resampled to a resolution of 0.87 x 0.87 x 0.87 mm (only changing in the craniocaudal direction), but the number of elements of the so generated FE model was too high to get reasonable computation times. Then, the images were resampled to a new resolution of 2.1 x 2.1 x 2.1 mm, using a linear interpolation of the grey level. In figure 5.2, the same slice of figure 5.1 is shown but now with the new resolution. As can be noticed the resolution is worse, mainly in the boundaries, but even the details are conserved, being this resample adequate given that a smoothing operation would be carried out afterwards to improve the FE model construction 105 Figure 5.2: Resampled image of the CT (2.1 x 2.1 mm resolution in the transverse plane). surfaces. As only one breast was modelled (the right one in this case) and no other parts of the body, the images were cropped in order to reduce the size of the model, and as a consequence, the number of elements, which is key in the computation time. Superiorly, the model was trimmed approximately at the height of the clavicle, and inferiorly below the inframammary fold. Medially, it was cut in the middle of the sternum (sagittal plane) and laterally by the coronal plane (see figure 5.3). The segmentation process started with an automatic segmentation with grey level thresholds. Then, a manual segmentation was performed to improve the automatic one and reduce errors. All the internal organs were removed (e.g. the heart or the lungs) and only the bones (ribs and sternum), muscle, skin and mammary tissue were left. Finally, a smoothing process with a recursive gaussian filter was carried out to smooth the boundaries. In figure 5.4 the result of the segmentation is shown. The bones appear in red color, the muscles in green and the fat in pink. The area shown in blue color was the brighter area and corresponds to the mammary gland, which is intermingled with adipose tissue. It was not possible to differentiate the fat from the glandular tissue in the CT image, so that both tissues had to be treated jointly. The properties of this mixture were calculated with the homogenized model presented in section 4.2, considering a 50% of fat, as estimated by the surgeons. It was not possible to differentiate the skin either, so it will be created later as a membrane covering the breast. This segmentation was then exported to the ScanFE module. The bones were assumed rigid and not included in the FE model. They were only segmented to locate the area in contact with the muscle, which was fixed in the FE model. Linear tetrahedral elements were selected to mesh the model and after meshing an additional smoothing was performed on the surface. Finally, the program generated an Abaqus input file. The skin was created by superimposing membrane elements on the exterior face of the adipose tissue 112 Finite element models of the breast tried to analyse their influence. In the model, the muscle was only in contact with the adipose tissue, and the modelling of this interaction was explained above. 5.4.3 Glandular tissue The glandular tissue was considered as rigidly bonded to the fat over its contact surface. It has no contact with the rest of the tissues, except for a small area in which it is connected to the skin (near the nipple area) and here it was considered rigidly bonded too. In the FE model in which the glandular tissue and the fat were modelled jointly, the boundary conditions of this mixture were equivalent to those of the adipose tissue, explained above. In figure 5.11, the glandular tissue is not shown because it is inside the model, surrounded by fat and skin. Therefore, this tissue is not involved in the four cut boundaries. 5.4.4 Skin As said above, the skin is closely bonded to the adipose tissue and, thus, it was modelled as perfectly bonded. With respect to the four cut boundaries, the skin of the model is connected to the rest of the skin and therefore, some kind of elastic bound must exist. In the medial boundary (in the sternum), the skin can be considered fixed mainly because it is very close to the bone and its movement is limited, due to the symmetry with respect to the sagittal plane. In the other three boundaries, different conditions were tested: fixed, totally free and an intermediate situation with springs, so simulating the effect of the rest of the skin. The precise conditions will be specified in the corresponding section. 5.5 Material models used All the tissues were considered near-incompressible, hyperelastic and isotropic (except for sections 5.6.12 and 5.6.13.3, in which a fibred model was used). In table 5.1, a summary of the the material models used through this chapter is presented. A short name is given to each in order to simplify their reference through the subsequent sections. The original constants of some material models were divided by certain factors in this thesis, as will be seen later. This change is denoted with a subscript in the name of the model, indicating the factor which divides the original constants. Moreover, when the homogenized model presented in section 4.2 was used, a superscript indicates the fat percentage considered. For example, the H50 100 material model is the homogenized model with a 50% of fat, with the constants divided by a factor of 100. FE Models 113 Name Material model Taken from Strain energy function Tissue Original material constants BM Rubin- Bodner Barbarino et al. [43] Ψ = µ0 2q(eqg −1), g= 2 m1J−1−ln(J)+(1−ω)m2(¯ I1−3) Muscle µ0= 3.7 MPa, q= 25, m1= 595, m2= 2.16 ·10−3,ω= 0 Ohl Linear elastic Ohl et al. [91] Muscle E= 790 kPa, ν= 0.45 Sf5 terms polynomial Samani et al. [41] Ψ = C10(I1−3) + C01(I2−3) + C11(I1− 3)(I2−3) + C20(I1−3)2+C02(I2−3)2Fat C10 = 310, C01 = 300, C11 = 2250, C20 = 3800, C02 = 4720 (Pa) T T First order Ogden Model proposed in this thesis, section 4.3 Ψ = µ α(λα 1+λα 2+λα 3−3) Fat µ= 67.869 (Pa), α= 8.575 Hp f 5 terms polynomial Homogenized model proposed in section 4.2 Ψ = C10(I1−3) + C01(I2−3) + C11(I1− 3)(I2−3) + C20(I1−3)2+C02(I2−3)2Mixture fat-gland Function of p(see section 4.2) BS Rubin- Bodner Barbarino et al. [43] Ψ = µ0 2q(eqg −1), g= 2 m1J−1−ln(J)+(1−ω)m2(¯ I1−3) Skin µ0= 1.7 MPa, q= 36, m1= 1294, m2= 8 ·10−4,ω= 0.9 P P 5 terms polynomial P´erez del Palomar et al. [18] Ψ = C10(I1−3) + C01(I2−3) + C11(I1− 3)(I2−3) + C20(I1−3)2+C02(I2−3)2Skin C10 = 31, C01 = 30, C11 = 22.5, C20 = 50, C02 = 60 (Pa) GOH Gasser- Ogden- Holzapfel Gasser et al. [46] Ψ = C10 2(I1−3) + k1 2k2ek2[κ¯ I1+(1−3κ)¯ I4−1]2−1Mixture fat-gland C10 = 300 Pa, k1= 1700 Pa, k2= 15.4, κ= 0.13 Table 5.1: Summary of the material models used in chapter 5. The subscript fstands for the factor dividing the original constants. The superscript pstands for the percentage of fat considered. 5.6 FE Models In this section, the different FE models which have been built are presented. Their advantages or problems are discussed and the subsequent modification introduced to each model will be justified in terms of improvement of the deformed shape compared with the 3D scan image. Only those models which are useful to understand the process followed during the development of this thesis are presented. Many other models, checks, attempts and convergence tests carried out are not included. In tables 5.2, 5.3 and 5.4, a summary of the FE simulations presented in this work is included. 5.6.1 Muscle fixed to the ribs and sternum, without skin. The first model was made with the simplest configuration. Thus, like in other breast models in the literature, all the materials were considered bonded to the surrounding tissues in the contact between each other, and their properties were taken from the literature. Only those nodes of the muscle in contact with the ribs and sternum were fixed, the rest were considered free to move. The skin was not modelled in these first models to analyse separately the behaviour of the other tissues. For the muscle, the BM model, presented in table 5.1, was used. The values of the constants are also shown in table 5.1. A UHYPER subroutine was programmed to implement the Rubin and Bodner model, as it is not included as a material model in Abaqus. The properties of the muscle are not the most influential ones, due to the fact that it is much stiffer than the breast tissues and that it is not much deformed with the position change. For the adipose tissue, the S1model was used, i.e., with the original properties, shown 114 Finite element models of the breast Section Skin Breast tissues Muscle Boundary conditions Other 5.6.1 No Fat (S1) and fat-gland mixture (Hp 1)BM Contact muscle-ribs/sternum fixed Supine-prone simulation, with 3 proportions for the fat-gland mixture (30%, 40% and 50% of fat) Supine-standing up simulation, with 3 proportions for the fat-gland mixture (30%, 40% and 50% of fat) 5.6.2 BS. Membrane Fat (S1) and fat-gland mixture (H50 1)BM Contact muscle-ribs/sternum fixed 5.6.3 BS. Membrane Fat (S1) and fat-gland mixture (H50 1)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed 5.6.4 BS. Membrane Fat (Sj) and fat-gland mixture (H50 j)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries 2 simulations using S100 -H50 100 and S200 -H50 200 BS. 3D solid Fat (S200) and fat-gland mixture (H50 200)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries No Fat (S200) and fat-gland mixture (H50 200)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries 5.6.5 No Fat (Sj) and fat-gland mixture (H50 j)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries 5 simulations using S200 -H50 200,S400 -H50 400,S800 - H50 800,S1600 -H50 1600 and S4000 -H50 4000 5.6.6 BS. Membrane Fat (Sj) and fat-gland mixture (H50 j)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries 2 simulations using S1-H50 1and S100 -H50 100 Fat (T T ) BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Fat fixed in medial and lateral boundaries 5.6.7 BS. Membrane Fat (S1) and fat-gland mixture (H50 1)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed Thin layer of tissue between muscle and breast tissues (neo-Hookean). 3 simulations varying C1 5.6.8 BS. Membrane Fat (Sj) and fat-gland mixture (H50 j)No Perimeter of skin fixed 4 simulations using S10 -H50 10 ,S100 -H50 100,S400 - H50 400 and S2000 -H50 2000 5.6.9 BS. Membrane Fat (S100) and fat-gland mixture (H50 100)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle and skin fixed. Spring connection muscle-breast tissues Several spring stiffness and distributions 5.6.10 PP. Membrane Fat (S100) and fat-gland mixture (H50 100)No Perimeter of skin fixed. Fat fixed in medial boundary BS. Membrane Fat (S100) and fat-gland mixture (H50 100)No Perimeter of skin fixed. Fat fixed in medial boundary 3 simulations changing µ0and 6 varying q Lateral and medial boundaries of skin fixed. In upper and lower boundaries, displacements constrained in craniocaudal direction. Fat fixed in medial boundary 3 simulations varying q, 4 varying m1and 3 changing m2 Table 5.2: Summary of the FE simulations presented in this work, ordered by section. FE Models 115 Section Skin Breast tissues Muscle Boundary conditions Other 5.6.11.1 PP. Membrane Unique, fat-gland mixture (H80 j)Ohl Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in medial-lateral and anteroposterior directions allowed, but not in craniocaudal 3 simulations using H80 10 ,H80 100 and H80 1000 Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum 2 simulations using H80 100 and H80 1000 5.6.11.2 PP. Membrane Unique, fat-gland mixture (H80 1000)Ohl Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In the lower boundary of skin, displacements constrained in craniocaudal direction. Springs in the lateral boundary of skin Simulations with several spring stiffness Contact in the whole surface between muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In the lower boundary of skin, displacements constrained in craniocaudal direction. Springs in the lateral and upper boundaries of skin Simulations with several spring stiffness 5.6.11.3 PP. Membrane Unique, fat-gland mixture (H80 1)Ohl Contact surface muscle-breast tissues in a stripe of the lateral boundary. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained 5.6.11.4 PP. Membrane Unique, fat-gland mixture (H80 j)Ohl Contact surface muscle-breast tissues in a stripe of the lateral, upper and lower boundaries. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained 3 simulations using H80 1,H80 10 and H80 100 Unique, fat-gland mixture (H80 100)Ohl Contact surface muscle-breast tissues in a stripe of the lateral, upper and lower boundaries. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained. Springs in the lateral and upper boundaries of skin Simulations with several spring stiffness Contact surface muscle-breast tissues in a stripe of the lateral, upper and lower boundaries. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained. Springs in the lateral and upper boundaries of skin, and in the whole contact muscle-breast tissues Simulations with several spring stiffness Table 5.3: Summary of the FE simulations presented in this work, ordered by section. 116 Finite element models of the breast Section Skin Breast tissues Muscle Boundary conditions Other 5.6.12 PP. Membrane Unique. GOH Ohl Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum Fibres: anteroposterior direction, with dispersion. 2 simulations changing C10 Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained Fibres: anteroposterior direction, without dispersion. 5 simulations changing C10 Contact in the whole surface muscle-breast tissues. Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in sternum. In upper and lower boundaries of skin, displacements in craniocaudal direction constrained. Springs in the lateral boundary of skin Fibres: anteroposterior direction, without dispersion. Simulations with several spring stiffness 5.6.13.1 BS. Membrane Fat (S1600) and fat-gland mixture (H50 1600)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in medial and lateral boundaries 2 simulations with different skin prestress Fat (S4000) and fat-gland mixture (H50 4000)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in medial, lateral and upper boundaries 2 simulations with different skin prestress Fat (S100) and fat-gland mixture (H50 100)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle and skin fixed. Breast tissues fixed in medial boundary. In upper and lower boundaries, displacements in craniocaudal direction constrained. No connection muscle-breast tissues Prestress only in half of the skin. 10 simulations with different skin prestress 5.6.13.2 BS. Membrane Unique, fat-gland mixture (H80 1)BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in medial boundary 4 simulations using H80 1,H80 100, H80 500 and H80 1000 Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in medial boundary. In upper and lower boundaries, displacements in craniocaudal direction constrained. No connection muscle-breast tissues 4 simulations using H80 1,H80 5, H80 10 and H80 100 5.6.13.3 BS. Membrane Unique. GOH BM Contact muscle-ribs/sternum fixed. Perimeter of muscle fixed. Breast tissues fixed in medial and lateral boundaries. In upper and lower boundaries, displacements in craniocaudal direction constrained Fibres aligned with the maximum principal Cauchy stress direction Table 5.4: Summary of the FE simulations presented in this work, ordered by section. FE Models 117 Figure 5.12: Deformed shape of the model for a 50% fat proportion superimposed with the real deformation in prone position, presented in section 5.6.1. in table 5.1. For the mixture of gland and fat, the H30 1,H40 1and H50 1models were used (see table 5.1). As estimated by the surgeons who treated the patient under study, her breast had about 50% of fat. Then 3 different fat percentages were tried for the mixture: 30%, 40% and 50%. The constants of the 5 terms polynomial model for each fat proportion are shown in table 5.5. Fat proportion C10 C01 C11 C20 C02 (Pa) 30% 324.53 288.26 3785.07 6468.84 7905.41 40% 322.38 290.71 3574.00 6104.66 7404.93 50% 320.10 293.31 3334.66 5656.75 6921.28 Table 5.5: Values for the mixture of fat and gland, obtained from section 4.2. Two static analyses were run using Abaqus Standard to simulate the change from supine to prone and standing up positions. The only applied load was the gravity. As the initial model was built in supine position, the gravity was applied twice in posteroanterior direction to obtain the prone deformed shape. In the first step to “eliminate” the gravity1and in the second step to simulate it in prone position. In the case of the standing up position, the gravity was first applied in posteroanterior direction to estimate the unloaded configuration, and then applied again in craniocaudal direction to obtain the standing up deformed shape. The accuracy of the deformation was evaluated quantitatively as explained in section 5.3. In figure 5.12 the deformed shape in prone position is shown in red superimposed with the real deformation (shown in green in the figure). Table 5.6 shows the mean absolute and mean relative errors for each fat 1The deformed shape of this first step could be thus considered a reference configuration (unloaded). This is not exactly true because this configuration is actually loaded (gravity in posteroanterior direction) and therefore it is not stress free. The problem comes from the fact that the model was built in supine position and therefore there were no loads or stresses in the model in this configuration, but in reality, this is a loaded state. This will be discussed later on (see section 5.6.13.2). 118 Finite element models of the breast 30% fat 40% fat 50% fat Mean absolute error 23.8±13.9 mm 23.6±14.4 mm 24.8±15.3 mm Mean relative error 86.2% 85.5% 89.8% Table 5.6: Quantitative comparison between the deformed shape of the model (presented in section 5.6.1) and the real deformation in prone position. Figure 5.13: Deformed shape of the model for a 50% fat proportion superimposed with the real deformation in standing up position, presented in section 5.6.1. proportion. Figure 5.13 presents the deformed shape in standing up position (in red) together with the real deformation (in green). In table 5.7 the mean absolute and mean relative errors are shown for each fat proportion. The errors were computed as explained before and depicted in figure 5.10. As can be noticed in figure 5.12, the simulated deformation was very far from the real deformation. To be precise, the model was almost not deformed, indicating that the behaviour was much stiffer than in reality. The errors in table 5.6 clearly show that the model behaviour was not correct in this case. Moreover the fat proportion had little influence, in the specified range and with the Hp 1models. In figure 5.13, corresponding to the standing up position, the simulated deformation was very far from the real one, and the model was almost not deformed, as in the previous simulation. However, it seemed that the deformed shape of the model was closer to reality than in the previous simulation. In view of table 5.7, it can be seen that the mean absolute error was lower than in the prone case, but only because the 30% fat 40% fat 50% fat Mean absolute error 5.8±3.1 mm 5.4±3.1 mm 5.3±3.1 mm Mean relative error 90.6% 84.4% 82.8% Table 5.7: Quantitative comparison between the deformed shape of the model (presented in section 5.6.1) and the real deformation in standing up position. FE Models 119 displacements were smaller. If the relative errors are compared, the percentages are of the same order of magnitude. The standing up simulation was carried out to check if the unrealistic deformed shape obtained for the prone position was related only with the change of position from supine to prone, or on the contrary, it was an inherent problem associated with the model itself. The results indicated that the simulated position was not the problem and that the discrepancies may arise from different reasons, but the model was much stiffer than the real breast and one of these reasons could be that the original Samani’s properties, used for the fat and from which the Hp 1models were obtained, were inadequate. Given that the relative errors were similar in both positions and it was easier to check visually whether the prone position was correct, in the subsequent simulations this was the simulated position and the visual check was the preferred option to evaluate the accuracy of the simulated deformed shape. 5.6.2 Muscle fixed to the ribs and sternum, with skin The model presented in the previous section (with the H50 1model for the fat-gland mixture) was run without the skin to analyse the behaviour of the other tissues. But, obviously, the skin should be included in the model as it is in the reality. The inclusion of the skin is the only difference with the model of the previous section. The skin was generated by selecting the outer faces of the tetrahedral elements and superimposing membrane elements over those faces such that both types of elements shared the outer nodes. The BS model was used for the skin with the material constants shown in table 5.1. As can be seen in the deformed shape of figure 5.14, the deformation was very small. The displacements were measured and compared with the previous model and were even smaller, as expected, since the addition of the skin made the model stiffer, obviously. Nonetheless, the skin was kept in the model from this point on (except for section 5.6.5). 5.6.3 All the boundaries of the muscle fixed, with skin In the previous models, even though the displacements were small and the global behaviour very stiff, the muscle was deformed in a non anatomical way in certain areas of the boundaries, being quite separated from its true position. For example, in the medial-inferior area, corresponding to the part below the sternum, as shown in figure 5.15. In the lateral boundary of the model, the same occurred. This strange deformations may be due to the lack of restrictions in the boundaries of the muscle, which were considered free to move, while in reality the pectoral muscles are joined to other big groups of muscles creating a strong bond. To avoid this separation of the muscle from its true position (e.g., from the sternum, as seen in figure 5.15), the muscle was also considered fixed in the whole perimeter of the model. The rest of parameters remained unchanged from the previous sections. The deformed shape is compared with that of section 5.6.2 in figure 5.16. It can be seen that they are almost identical except for the area close to the sternum and the 120 Finite element models of the breast Figure 5.14: Deformed shape of the model with skin in prone position, presented in section 5.6.2. Figure 5.15: Deformed and undeformed shapes of the model with skin in prone position, presented in section 5.6.2. Lateral view of the sternum. FE Models 121 Figure 5.16: Left: model with the muscle fixed only in the contact with ribs and sternum, presented in section 5.6.2. Right: model with the muscle fixed in the whole perimeter of the model, presented in section 5.6.3. Upper view. lateral boundary. These strange deformations of the muscle were eliminated but, as expected, the model was slightly stiffer than the previous ones, because of the new boundary conditions imposed. 5.6.4 Sensitivity of the model to the stiffness of the fat-gland properties, with skin In the true deformation of the breast, the soft tissues are displaced towards the axilla in supine position, whereas in prone position they move closer to the sternum. After realizing in the previous models that they were much stiffer than the reality, even without skin and with different boundary conditions, a sensitivity analysis was performed on the elastic properties assumed for the fat and glandular tissues. For this purpose, the material parameters of the strain energy function were divided by certain factors to analyse their effect on the deformed shape. Apart from dividing the material constants by a certain factor, some new boundary conditions were added to this set of simulations, to bring the behaviour of the model closer to reality. The adipose tissue was fixed in the medial (sternum) and lateral boundaries. In the sternum, it was fixed because in previous simulations, as well as in reality, it was checked that the displacement of this area was negligible. Regarding the lateral boundary, it was observed that if the adipose tissue was free, the skin of this boundary dropped excessively and rotated in an incorrect direction, while in the true deformation the skin remained near the chest, as can be noticed in figure 5.8. Probably, this deformation would be constrained by the rest of the skin, not included in the model. For the moment, the traction that the skin outside the model is exerting 128 Finite element models of the breast Pa, some convergence problems arose1. The layer of material was very thin and when large deformations were produced, some elements became very distorted. This model simply represent a transition to the following models, but its performance was the same as in the previous simulations, because there was nothing in the model that allowed the lateral movement of the breast. 5.6.8 Complete lack of connection between the muscle and the breast tissues The connection between the muscle and the breast tissues can play an important role in the behaviour of the breast. In the previous models both tissues were rigidly bonded. In this section the other extreme case was studied: muscle and breast tissues were not bonded and no restriction or connection existed between both tissues. This approach has a severe limitation: a separation of the surfaces of both tissues, which were initially in contact, would create a void which would be hard to explain from a physiological point of view. Therefore, it only makes sense as another extreme situation in which the retromammary space was filled with a material of null (or almost zero) stiffness. As the change of position from supine to prone was simulated, the role of the muscles had no influence on the results and they were eliminated from the model. The boundaries of the skin were then fixed, to support the weight of the breast tissues. The BS material model was used to characterise its behaviour. Four models with different material constants for the adipose tissue and the fat-gland mixture were used to study the evolution of the deformed shape: S10-H50 10 ,S100-H50 100,S400-H50 400 and S2000-H50 2000. The deformed shape of the model using S10-H50 10 is shown in figure 5.23 (right), next to the 3D scan image (left). They are not superimposed for clarity. In figure 5.24 the models using S100-H50 100 and S2000-H50 2000 are shown together. The model using S400 for the fat and H50 400 for the fat-gland mixture is not shown, because it followed the same tendency as the shown models. As can be seen in figure 5.23, the displacements of the model using the S10 material model for the adipose tissue and the H50 10 for the fat-gland mixture were larger than in the reality. In fact, the breast behaved similarly as a balloon filled with a fluid, hanging from the boundaries and forming only one bulge. In the real breast, the tissues are displaced near the sternum, forming one bulge, wide in craniocaudal direction and narrow in mediolateral direction. Moreover, the skin is not separated from the chest until the bulge is formed near the sternum, causing two pronounced changes in the curvature, responsible for its particular shape. These differences with reality can be observed in figure 5.23. In the deformed shape of the model, these two areas with different curvatures did not appear, because the breast tissues were completely separated from the muscle. Moreover, the lowest point of the deformed shape was far from the sternum, because there 1To try to overcome the difficulties generated by the distortion of the elements, the capabilities of the adaptive mesh in Abaqus Standard were explored. However, it is not possible to use adaptive mesh with tetrahedral elements in Abaqus Standard. Besides, it is indicated for wear and ablation processes, not for large deformation ones. In Abaqus Explicit it is not possible to use the adaptive mesh either. All the tetrahedral and hexahedral elements are full integrated or reduced integrated with hourglass control, and adaptive mesh does not support full integration or hourglass control elements. Then, the adaptive mesh procedure was abandoned. FE Models 129 Figure 5.23: Model without connection between the muscle and breast tissues. Left: 3D scan of the breast in prone position, highlighting the differences between the real breast and the FE simulated deformed shape. Right: model using S10 for the adipose tissue and H50 10 for the fat-gland mixture. Upper view. Section 5.6.8. Figure 5.24: Model without connection between the muscle and breast tissues. Left: S100 and H50 100. Right: S2000 and H50 2000. Upper view. Section 5.6.8. 130 Finite element models of the breast Figure 5.25: Left: 3D scan of the breast in prone position. Right: model with springs connecting the muscle with the breast tissues. Upper view. Section 5.6.9. was nothing in the breast model that could cause the breast to come closer to the sternum. In figure 5.24, it is clearly seen that the deformation and displacements in the right-hand side model were really enormous, even compared with the left-hand side model, which also presented larger displacements than the reality. The simulation with this so flexible material was done to analyze the evolution of the deformed shape as the stiffness was lowered. In view of these results, it was concluded that the boundary conditions cannot be so flexible, and that there must be something in the breast which avoids these large deformations. In section 5.6.6, it was seen that the S100 and H50 100 material models should be used in the FE breast model to obtain a similar behaviour to that using the material model experimentally obtained in this thesis. With the boundary conditions simulated here, using S10 and H50 10 (stiffer than S100 and H50 100) the deformations in the model were larger than those in the reality. Therefore, it seemed that the issue was more related to the boundary conditions than to the material properties of the tissues. 5.6.9 Spring connection between the muscle and the breast tissues In the previous section, the displacements of the breast model were larger than those of the real organ for the reference material constants (S100 and H50 100). Next, a model with a spring connection between the muscle and breast tissues was built to reproduce an intermediate situation and so to mimic a structure of loose connective tissue or weak ligaments. Therefore, in this model the muscle is included again using the BM material model, fixing the nodes in contact with the ribs and those of the four cut boundaries. The muscle and breast tissues were not bonded, but each node of the muscle surface was connected by a linear spring with the facing node of the posterior breast surface (one to one correspondence), producing a force FE Models 131 proportional to the distance between both nodes 1. The skin was modelled with the BS material and fixed in the four boundaries again. For the fat, the S100 material model was used and for the mixture of fat and gland, the H50 100. The stiffness of the springs was set equal to 0.1 N/m, with a density of 7.5 springs/cm2, which results in an equivalent stiffness of 7.5·10−3N/cm/cm2. This value is not based on the literature, but taken to obtain approximately the same displacement in anteroposterior direction as this of the 3D scan image of the real breast in prone position. The resulting deformed shape is shown in figure 5.25. Again, it was not satisfactory. Another models with springs were tried, dividing the whole surface in four areas (upper-medial, upperlateral, lower-medial and lower-lateral) of approximately equal size. A different stiffness was assigned to the springs of each area. Different combinations of stiffness were tried (ranging from 0.01 to 1 N/m). The objective was to simulate a situation in which the stiffness of the retrommary space was different, depending on the region. The results improved, but they were not satisfactory yet and the number of variables that should be adjusted were increased. Obviously, dividing the surface in enough zones with springs with different stiffness, and even also dividing the breast tissues in areas with different properties, it would be possible to reproduce the shape of the breast. Nevertheless, the number of variables needed to be adjusted would be very large and the process would completely lose the physical sense. The breast behaviour would be totally faked. 5.6.10 Analysis of the influence of the skin properties In this section, the influence of the skin properties in the behaviour of the breast model is studied. The FE model without any connection between the muscle and the breast was used, because the properties of the skin played a crucial role in it. Therefore, the muscle was not included. For the breast tissues, the model S100 for the adipose tissue and the H50 100 for the mixture of fat and gland were used. They were fixed only in the medial boundary. 5.6.10.1 Comparison of two material models of the skin First, two different material models of the skin were compared: the BS model used in the previous simulations, and the PP model used by several authors [18, 21, 23] to cover the breast. The skin was fixed in its four boundaries. Figure 5.26 compares the deformed shapes of both models. They were practically equal (being the PP model slightly stiffer), indicating that both models behave similarly under the simulated conditions. Therefore, the BS model was selected to complete the rest of the study. 1Opposite to other options in which the spring force is proportional to the distance in a certain direction 132 Finite element models of the breast Figure 5.26: Model without connection between the muscle and breast tissues and material models S100 and H50 100. Left: Rubin-Bodner model for the skin. Right: polynomial model for the skin. Upper view. Section 5.6.10. 5.6.10.2 Sensitivity analysis of the Rubin and Bodner’s material constants A sensitivity analysis of the FE model to the values of the material constants of the BS material model was carried out. This model, as seen in section 3.2.5.1, depends on five material constants. The original values provided by Barbarino et al. [43] for the skin are included in table 5.1. Each constant was varied at a time ceteris paribus. The sensitivity analysis was only conducted for the first four constants. The last one, ω, is age related and was maintained with its original value. Moreover, as it multiplies m2in equation (3.31), both variables can be considered together for the purpose of this analysis. Analysis of µ0 The values tested for µ0were: 1.7 (initial), 17 and 85 MPa. In these simulations, the four boundaries of the skin were fixed. In figure 5.27 the deformed shapes of the three models with different values of µ0are compared. As can be seen, the model became stiffer as the constant µ0increased, reducing the maximum displacement, though with no effect on the global shape of the breast since all the nodes were proportionally displaced. Moreover, it can be seen that µ0needs to be considerably varied to notice changes in the deformed shape. Analysis of q Here, two boundary conditions were changed with respect to the previous simulations. In the real breast (figure 5.9, lateral and front views) the tissue moves towards the sternum forming a bulge which extends along the sternum in craniocaudal direction, while in the model there was a minimum in the center of the FE Models 133 Figure 5.27: Left: Original properties of the skin (µ0= 1.7). Centre: µ0= 17. Right: µ0= 85. Upper view. Section 5.6.10. Figure 5.28: Front view of the model without connection between the muscle and the breast tissues. Section 5.6.10. space delimited by the four line boundaries (as shown in figure 5.28). Therefore, in the new model, the upper (clavicle) and lower (closer to the abdomen) boundary conditions of the skin were changed: the displacements in the mediolateral and the anteroposterior directions were allowed, but not in the craniocaudal direction. That is to say, the nodes at the boundaries can only move within the plane of figure 5.27. The objective was to test the behaviour of the skin with a less constrained movement in these directions, allowing the displacement of the breast tissues towards the sternum. First the simulations were run with the same boundary conditions as in the previous section (analysis of µ0), i.e., with the four skin boundaries fixed. The values tested for qwere: 36 (initial), 360, 3600, 36000, 360000 and 600000. Then, both boundary conditions explained in the above paragraph, upper and lower, were changed for three values of q: 36, 36000 and 600000. The effect of the change of boundary conditions was also analysed. 134 Finite element models of the breast Figure 5.29: Four boundaries of skin fixed. Left: Original properties of the skin (q= 36). Centre: q= 36000. Right: q= 600000. Upper view. Section 5.6.10. Figure 5.30: Comparison of boundary conditions for the original properties of the skin. Left: four boundaries fixed. Right: upper and lower boundaries changed. Upper view. Section 5.6.10. Figure 5.31: Upper and lower boundaries changed. Left: Original properties of the skin (q= 36). Centre: q= 36000. Right: q= 600000. Upper view. Section 5.6.10. FE Models 135 In figure 5.29, the deformed shapes of three models (without changing the boundary conditions with respect to the previous simulations) with different values of qare compared: 36 (original value), 36000 and 600000. The rest of the models are not shown because they followed the same tendency and did not provide additional information to the deformed shapes shown in figure 5.29. The model was stiffer and the displacements smaller as the constant qincreased, and the deformed shape had its lowest point located more or less in the center of the model and without the pronounced curvatures seen in the real breast. In figure 5.30, both models for the original properties of the skin but changing the boundary conditions are compared. As can be seen, the effect on the deformed shape in craniocaudal direction was very important, because the upper and lower boundaries of the skin were released in mediolateral and anteroposterior direction in the model shown on the right-hand side. Seeing the model from an upper view, with these new boundary conditions, the shape of the breast was more similar to the reality in craniocaudal direction. Nonetheless, if the model is seen from a lateral view, both upper (sternum) and lower (abdominal) boundaries became excessively separated from the chest wall. Moreover, the displacement of the model was slightly larger in anteroposterior direction, but the lowest point remained far from the sternum (approximately at the same distance that in the model with the four boundaries of the skin fixed) and the pronounced curvatures were not produced. In figure 5.31 the deformed shapes (changing both new boundary conditions) for three values of qare compared: 36 (original value), 36000 and 600000. Similarly to figure 5.29, the model was stiffer as the constant qincreased. Comparing the models with different boundary conditions for the same value of q (figures 5.29 and 5.31), the maximum displacements in anteroposterior direction were slightly higher when the upper and lower boundaries were released. Analysis of m1 The values tested for m1were: 1294 (initial), 129.4, 12940 and 129400. The effect of varying the parameter m1is shown in figure 5.32 (the case m1= 12940 is not shown, as it did not provide additional information to the other images). Comparing the simulation for m1= 129.4 and m1= 129400, the deformed shape was approximately the same. The variation of this parameter within that range did not affect much the global behaviour. In the strain energy function (recall equations (3.30) and (3.31)), this parameter multiplies the terms which forces the incompressibility. Thus, as the skin is being modelled as incompressible, this term should be approximately zero and, therefore, the variations of m1has a very small effect on the behaviour. Analysis of m2 Finally the values tested for m2were: 0.0008 (initial), 0.00008 and 0.008. The effect of varying the 136 Finite element models of the breast Figure 5.32: Left: Original properties of the skin (m1= 1294). Centre: m1= 129.4. Right: m1= 129400. Upper view. Section 5.6.10. parameter m2is shown in figure 5.33. As can be noticed, the deformation did not vary much. Similarly to the case of m1, in the strain energy function (recall equations (3.30) and (3.31)) m2multiplies the term (1 −ω) being ω= 0.9 as proposed by Barbarino et al. [43]. Moreover, the original value of m2is very small. Because of this two reasons, the effect of m2on the behaviour is noticeable but low. 5.6.11 Contact surface between the muscle and the breast tissues The previous models presented in this chapter failed mimicking the real behaviour of the breast. On the one hand, in the FE models with a rigid bond between the muscle and breast tissues, the latter were not displaced to the sternum and the deformations using the material properties provided by the literature were very small. Making the tissues more flexible did not solve the problem and even created a new one, because two bulges appeared in the model. On the other hand, when the muscle and the breast tissues were not bonded, the deformed shape was more similar to the real one, but the tissues were far from the sternum and the two areas with different curvatures (see figure 5.23) could not be distinguished. Based on the anatomy, the experimental observation and all the simulations carried out, the deformation of the real breast in prone position seemed to be caused by a sliding of the breast tissues over the muscle. The retromammary space located between both tissues is formed by loose connective tissue, ligaments and a fascial structure. Therefore, it is possible that this anatomical structure allows a finite sliding between the two types of tissues over their contact surface. For this reason, a model with a contact interaction between the muscle and the breast tissues was built. The contact was modelled with a finite sliding contact formulation, surface to surface, without friction and with no separation, forcing both surfaces to remain in contact. A different FE model to the one used in the FE Models 137 Figure 5.33: Left: Original properties of the skin (m2= 0.0008). Centre: m2= 0.00008. Right: m2= 0.008. Upper view. Section 5.6.10. previous sections was used for these simulations, the last one presented in section 5.2, with 138229 elements. It is important to note that this model is the first one in which the adipose tissue and the fat-gland mixture were considered as a unique material (a total fat-glad mixture), modelled with the H80 jmaterial model. This model was constructed to smooth the surface between the muscle and the breast tissue, because the previous model did not converge under gravity loads (even in incremental load analyses for low percentages of the applied gravity), as the surface was not smooth enough. Three different types of contact were simulated: 1) with the whole surface between muscle and fat-gland mixture in contact, 2) with only the lateral line (the boundary line) of the fat-gland mixture enforced to glide on the muscle and 3) with the lateral, upper, and lower boundary lines of the fat-gland mixture enforced to glide on the muscle. Therefore, in the last two contact interactions, the rest of the surface of the fat-gland mixture can lose contact from the muscle, except for the lines in which the no separation contact is modelled. The simulations were carried out with Abaqus Explicit, because the surfaces did not present a regular shape even after smoothing, thus making difficult the convergence with Abaqus Standard. With Standard, the simulations neither were successfully finished nor were aborted with a substantial percentage of the gravity load applied, even after trying different combinations of parameters of the contact formulation. In Abaqus Explicit only dynamic simulations are allowed. As the deformed shape under static conditions was desired, some simplifications were implemented to obtain it: the load was applied slowly, damping was added in order to reduce the oscillations around the static deformation and the simulation time was long enough to allow time for the model to tend to the static deformed situation. In the previous models used in Abaqus Standard, the skin and the muscle tissues were modelled with a Rubin and Bodner strain energy function by means of a user subroutine (UHYPER). But, in Abaqus Explicit, 144 Finite element models of the breast Figure 5.40: Model with contact surface between the muscle and breast tissue, this last one modelled with the GOH model, with fibres in anteroposterior direction and with dispersion. Left: 3D scan in prone position. Centre: C10 = 300 Pa. Right: C10 = 1 Pa. Upper view. Section 5.6.12. function of the modified invariant ¯ I4, represents the anisotropic contribution of the fibres, modelling the strong stiffening of the collagen fibres when stretched. The preferential direction of the family of fibres is defined by the unit vectors a0, with respect to the material configuration. The material constants of the model were taken from [47] and shown in table 5.1 as GOH model. Two cases were studied in Abaqus Explicit: •With fibre dispersion (κ= 0.13): the parameter C10 is related to the stiffness of the matrix where the fibers are embedded. The values 300 (original) and 1 Pa were used for the constant. The muscle was fixed in its four cut boundaries and in the contact with the ribs. The breast tissue and the skin were fixed only in the sternum, being the rest of the boundaries free to move. The whole surface between the muscle and the breast tissue was modelled as a contact surface, not allowing the separation between both surfaces. •Without fibre dispersion (κ= 0): the parameter C10 was also varied: 600, 400, 200, 100 and 20 Pa. k1was changed with respecto to the original value, to allow a higher deformation (k1= 850 Pa). The boundary conditions are the same as in the previous case with fibre dispersion, except for the upper and lower boundaries (clavicle and abdomen) of the skin, which were not totally free to move in these simulations. More precisely, the movement in craniocaudal direction was constrained. In figure 5.40 two deformed shapes of the model with fibre dispersion, for C10 equal to 300 and 1 Pa, are shown and compared with the 3D scan image in prone position. In figure 5.41 the deformed shape of the model without fibre dispersion is shown, for the values of C10 600 and 200 Pa. In figure 5.43, the upper and front views of the deformed shape of the model without fibre dispersion is shown for C10 = 100 Pa. The models for the other values of C10 are not shown because they do not provide additional information to that presented in the previous figures. FE Models 145 Figure 5.41: Model with contact surface between the muscle and breast tissue, this last one modelled with the GOH model, with fibres in anteroposterior direction, without dispersion. Left: 3D scan in prone position. Centre: C10 = 600 Pa. Right: C10 = 200 Pa. Upper view. Section 5.6.12. Figure 5.42: Model with contact surface between the muscle and breast tissue. Left: isotropic model for the breast tissues. Right: anisotropic model for the breast tissues. Upper view. Section 5.6.12. As can be seen in figure 5.40, the deformed shape is quite unrealistic, specially in the case with C10 = 1Pa (right) and likely due to numerical problems. The same simulations were conducted but without fibre dispersion. The objective was to check if the resulting deformed shape is more realistic without fibre dispersion. The deformations were normal, as can be seen in figure 5.41, and thus the fibre dispersion was thought to be responsible for the numerical problems. In figure 5.42, the deformed shape on the right-hand side of figure 5.41 is compared with the deformed shape of the model with a contact surface between the muscle and the breast tissues, but with the latter modelled as isotropic (the same model presented in section 5.6.11.1 in figure 5.35). It can be seen that the displacement of the anisotropic model in anteroposterior direction is smaller, but in mediolateral direction is even larger. Besides, the lowest point of the breast was not as markedly abrupt as in the previous cases 146 Finite element models of the breast Figure 5.43: Model with contact surface between the muscle and breast tissue and C10 = 100 Pa. Left: upper view. Right: front view. Section 5.6.12. with an isotropic material. Therefore, the material was assumed slightly more flexible, to try to reproduce a closer deformation to the real one, more displaced to the sternum and with a slightly larger displacement in anteroposterior direction. For C10 = 100 Pa, the simulation converged but the deformed shape had many wrinkles as can be seen in figure 5.6.12, mainly in the front view. The simulation did not converge for C10 = 20 Pa. To solve the convergence issue, some solutions were tried: changing some parameters, as the ones which control the damping, and adding springs with different stiffness (0.05 and 0.1 N/m) in the boundary to stabilize the movement, but the convergence was not reached with these strategies. 5.6.13 Other models This section presents the results of some simulations carried out with other models. These models are not less important than the rest. The reason why they are presented apart is that they do not fit well in the time line which has been followed during this chapter and it was preferred to explain them separately. 5.6.13.1 Skin prestress The skin prestress was first studied because of two reasons. The first one is that it really exists in the skin of the body as will be discussed later on. Moreover, when convergence problems arose due to skin compression (section 5.6.4), the skin prestress was thought to play an important role in the simulations. The skin can experience unstable behaviours such as buckling under compression, in the body and particularly in the breast (for example if one pinches a piece of skin between the fingers), but at least in the organ studied FE Models 147 here, that is not a normal situation. However, a certain compression appeared in the simulations. This compression could actually be a reduction of its tension if the skin was prestressed enough. The second reason is that it was also known that the initial configuration of the model in supine position is not a stress free state and that, maybe, it had an important influence on the results of the FE model, but this last point will be explained better in section 5.6.13.2. Some information can be found in the literature about skin prestress in the human body, and particularly in the breast. However, the majority of the articles focus on Langer’s lines, their distribution and the magnitude of the skin retraction when cut [126, 127, 128, 129], but few give a quantitative value of the prestress. Moreover, the values differ quite much among the different authors. For example, Alexander and Cook [130] used a suction device to investigate the prestress of the skin in the back. They found that this prestress was between 5.6 and 24.8 kPa. Flynn et al. [131] did experimental tests in the forearm and fitted them with FE analyses, determining that the prestress was in the range between 28 and 92 kPa. In the same year, Flynn et al. [132] did more experiments with the same procedure, but giving a different range for the prestress, now between 1.1 and 7.6 kPa. Diridollou et al. [133] did also experimental tests in the skin of the forearm, estimating an average prestress of 13.5 kPa. Therefore, it is clear that the prestress is a real characteristic of the skin, but its real magnitude is not known with certainty. Uniform prestress As presented in section 5.6.4, the convergence problems due to compression in the skin appeared when a model with the following conditions was simulated: muscle fixed in the four cut boundaries and in its contact with the ribs, breast tissue and skin fixed in the medial (sternum) and lateral boundary, muscle and breast tissues rigidly bonded and using the Sfmaterial model for the fat and the H50 fmodel for the fat-gland mixture, with f > 100. For the skin and the muscle, the BS and BM models were used respectively. Lets us recall (see section 5.6.4) that the same model without skin converged without problems. For example, in figure 5.44 on the left-hand side, the model without skin and with the S1600 material model for the fat and the H50 1600 for the fat-gland mixture is shown again. For this reason, this model, now with skin, was used to study the skin prestress and the convergence. The prestresses were applied in an equibiaxial way, i.e., the same magnitude was applied in two perpendicular directions of the membrane elements of the skin, creating a uniform prestress. With the S1600 and H50 1600 material models, prestresses of 100 and 150 Pa were introduced in the model before the gravity loads were applied. With the S4000 and H50 4000, skin prestresses of 150 and 1250 Pa were studied. The boundary conditions of the skin were changed in these simulations, now fixing the upper and lower boundaries, to avoid the release of the pretension in craniocaudal direction. In figure 5.44, the deformed shape using the S1600 and H50 1600 material models and a prestress of 150 Pa 148 Finite element models of the breast Figure 5.44: Analysis of the effect of skin prestress. Left: S1600 and H50 1600, without skin. Centre: S1600 and H50 1600, with skin prestress (150 Pa). Right: S4000 and H50 4000, with skin prestress (1250 Pa). Section 5.6.13.1. is shown in the center. On the right-hand side, the model with S4000 and H50 4000 and a prestress of 1250 Pa is included. On the left-hand side, the deformed shape of the model without skin and with the S1600 and H50 1600 models presented in section 5.6.4 is added to compare the results with those simulated in this section. The model using S1600 and H50 1600 and a uniform prestress of 100 Pa did not converge but few steps remained to complete the simulation (90% of the gravity applied). The prestress was not enough to avoid the compression in the skin. With 150 Pa of prestress, the simulation converged. Comparing the deformed shape of this model (in the centre of figure 5.44) with the left one in the figure, it can be seen that the deformation is almost the same. Therefore, the model without skin, and the one with a small skin prestress behaved similarly. Moreover, the prestress was able to overcome the convergence problem. In most of the areas of the model with skin prestress the surface was smoother, and the lowest point of the deformed shape slightly shifted towards the sternum, though not much. Besides, the displacement in anteroposterior direction was slightly larger, which at first seemed to be strange. The reason was that, as the skin was not fixed in the upper and lower boundaries, the prestress in craniocaudal direction could be relaxed, concentrating the breast tissue near the center and therefore producing a lower peak of the deformed shape. This boundary condition was not changed at first to allow a direct comparison with the previous model (on the left), but it was changed in the next model. The constants were divided by 4000 in the next simulations to allow a higher displacement in anteroposterior direction. Using 150 Pa of prestress, the simulation did not converge (it reached approximately the 50% of the applied load). With 1250 Pa, the simulation converged. Due to the fact that the tissues were more flexible in these simulations, a higher prestress was needed. As can be seen in figure 5.44 on the right, the effect of the new boundary conditions is clear, maintaining fixed the skin in the upper and lower boundaries. As noticed, the displacements in anteroposterior direction were similar compared with the other two models. However, the surface was smoother and the deformation near the boundaries was more normal. FE Models 149 Different prestresses in two regions of the skin Additionally, another model was built differentiating two areas of the skin separated by a central line along the craniocaudal direction, i.e., a half near the sternum (medial) and a half near the axilla (lateral). The objective was to analyse the effect of having different prestresses in each area. In the real breast, in supine position (initial configuration) the skin would be stretched due to the weight of the breast tissues. Moreover, the tension in the skin would not be uniform, as each portion of skin is not bearing the same load. For example, the skin near the sternum would have a higher prestress, whereas the skin near the axilla would be almost unloaded. This effect was modelled assigning a different prestress in different areas of the skin. So, the lateral half was not prestressed, while the medial one was. In these simulations, the muscle was fixed in the four cut boundaries and in its contact with the ribs and the skin was also fixed in the four boundaries. For the breast, besides fixing it in the medial boundary, the displacements in craniocaudal direction were also constrained in the upper and lower boundaries. The breast tissues and the muscle were not bonded in these simulations. For the breast tissues properties, the S100 material model was used for the fat and the H50 100 for the fat-gland mixture. The simulations were done using prestresses from 0 to 15 kPa. In figure 5.45, three deformed shapes are shown, with different prestress in the medial and lateral halves, corresponding to prestresses of 0, 500 and 5000 Pa. The rest of the prestresses below 5000 Pa are not included because they did not provide additional information to the shown deformed shapes. It can be seen that the skin prestress had an influence on the maximum displacements, and that as it increased, the breast tissues were slightly displaced to the sternum, although it is difficult to notice this effect with the naked eye. It is important to note that the magnitude of the prestress was considerably high compared to the low effect which it caused in the lateral movement of the tissues. For higher values, 10000 and 15000 Pa, the simulations did not converge due to element distortion. In conclusion, the prestress solved the convergence problem and had an effect on the deformed shape, but not important enough to obtain a deformation similar to the real one. 5.6.13.2 Pullback algorithm The FE models presented in this thesis were built from CT images taken in supine position, which was considered the initial configuration. Many authors who work with FE biomechanical models consider the medical images as the initial configuration, that is, stress free, although obviously it is not. In the case of the breast in the supine position, stresses and strains exist because the gravity is acting in this position, and they may play and important role in the subsequent simulations. Therefore, it would be good to calculate a reference state for the breast in which no external loads are applied. In this particular case, the objective was to “eliminate” the gravity from the model in supine position. Starting from this reference configuration and applying the gravity in different directions, the prone, supine and standing up positions would be obtained. 150 Finite element models of the breast Figure 5.45: Analysis of the effects of different skin prestresses as a function of the region. Left: prestress = 0 Pa. Centre: prestress = 500 Pa. Right: prestress = 5000 Pa. Section 5.6.13.1. In the previous simulations carried out in this thesis, the prone position was obtained applying twice the gravity in posteroanterior direction. With the first gravity load applied, it was intended to obtain a reference configuration, but it was not stress free. The approach presented in this section is closer to reality. To calculate the reference configuration of the breast, the algorithm proposed by Riveros et al. [120] was used. This is an iterative algorithm in which the reference configuration is updated in each step. It stops when, after applying the exiting loads to the updated reference configuration, the geometry of the medical images is obtained. In this particular case, it stopped upon the application of the gravity loads to the reference configuration, when the deformed shape in supine position, that taken from the medical images, was obtained (of course, with a given tolerance). The steps of the algorithm in the case of the breast are explained, although the reader is referred to Riveros et al. [120] for a deeper insight in the algorithm. The position of the nodes of the FE model mesh is called Xin a general basis. Ximage represents the configuration of the medical images. The objective is to obtain Xref , a mesh such that when the gravity loads are applied, the configuration Ximage is obtained. The reference configuration is updated in each step, Xk ref . 1. The reference configuration (configuration “without gravity”) is initialized with the supine position configuration. X0 ref =Ximage . 2. Gravity is applied in anteroposterior direction to the reference configuration of the (k−1) iteration, Xk−1 ref , to obtain the corresponding deformed configuration Xk def . FE Models 151 Figure 5.46: Scheme of the gravity load application. 3. Once the simulation is finished, the error is compared with the tolerance in each node iof the FE mesh. |Ximage −Xk def |< tol . 4. If the previous condition is not accomplished for every node, the reference configuration is updated. Xk ref =Xk−1 ref + (Ximage −Xk def ) . 5. Steps 2, 3 and 4 are repeated until condition 3 is fulfilled for every node. 6. If condition 3 is fulfilled, Xk ref is taken as Xref . To implement this algorithm, a Fortran script was written. This script called Abaqus to run the simulations, extracted the results, checked the condition to finish the loop and made the calculations of step 4. Abaqus Standard was used in these simulations. Once Xref was obtained, gravity was applied in posteroanterior direction to simulate the prone position. First, the pullback algorithm was used in the model in which the muscle and the breast tissue were rigidly bonded. The muscle was fixed in the four cut boundaries and in the contact with the ribs, and the breast tissue and the skin were fixed in the sternum. For the skin, the BS model was used, and for the muscle, the BM model. For the adipose tissue, the S1,S100,S500 and S1000 material models were used and for the mixture of fat-gland, the H50 1,H50 100,H50 500 and H50 1000. Then, the algorithm was applied to the model with no connection between the muscle and the breast tissue. The muscle was fixed in its four cut boundaries and in the contact with the ribs. The skin and the breast tissue were fixed in the sternum, and their displacements in the upper and lower boundaries were constrained in the craniocaudal direction. For the skin, the BS model was used, and for the muscle, the BM model. The S1and H50 1models were used for the fat and fat-gland mixture respectively. The results of the pullback algorithm for the model with the muscle and breast tissues rigidly bonded and the original Samani’s properties are presented in figure 5.47. On the left-hand side, the reference configuration is shown, i.e., the configuration ”without gravity” which was obtained after applying the pullback algorithm. In the center, the prone position obtained starting from this reference configuration is attached. Finally on the right-hand side, the prone position obtained without the pullback algorithm (that is to say, applying 152 Finite element models of the breast Figure 5.47: Comparison of using or not using the pullback algorithm. Left: reference configuration (“without gravity”). Centre: deformed shape in prone position using the pullback algorithm. Right: idem not using the pullback algorithm. Upper view. Section 5.6.13.2. Figure 5.48: Comparison of using or not using the pullback algorithm, with S100-H50 100. Left: reference configuration (“without gravity”). Centre: deformed shape in prone position using the pullback algorithm. Right: idem without using the pullback algorithm. Upper view. Section 5.6.13.2. Figure 5.49: Comparison of using or not using the pullback algorithm, without connection between the muscle and the breast tissue. Left: reference configuration (“without gravity”). Centre: deformed shape in prone position using the pullback algorithm. Right: idem using the pullback algorithm. Upper view. Section 5.6.13.2. FE Models 153 twice the gravity to the initial model in supine position, as in the previous simulations in this thesis) is shown. The results for S100-H50 100 are presented in figure 5.48. The results for the model with no connection between the muscle and the breast tissues are shown in figure 5.49. As in the previous figures, the reference configuration is shown on the left-hand side. In the center, the prone position using the pullback algorithm is included, and on the right-hand side, the prone position obtained without the algorithm is presented. In figure 5.47, comparing both models in prone position, the displacements were larger in the model using the pullback algorithm, although, due to the high stiffness of the model, it was difficult to notice the difference with the naked eye. Comparing both models in prone position in figure 5.48, it is clearly seen that in the deformed shape obtained using the pullback algorithm the displacements were larger in general and in particular the lateral movement. That means that the pullback algorithm produces a deformed shape closer to reality, with the same material, properties and boundary conditions. The process was repeated with the constants divided by 500 and 1000, to check what happened with a more flexible model, but the algorithm did not converge due to bad shaped or even zero volume elements appearing while updating the reference configuration during the iterative process. Regarding the models without connection between the muscle and the breast tissues, compared in figure 5.49, it can be noticed that in the new model (center) the breast tissue shifted more in anteroposterior direction and to the sternum that in the model without the pullback algorithm, therefore improving the results. However, the lateral boundary became very separated from the muscle surface, as can be seen in figure 5.49, and this is not very realistic. Using a model with a contact surface between the breast tissue and the muscle together with the pullback algorithm was not possible. In the first step of the algorithm, both surfaces were in contact. However, when the reference configuration was updated, the surfaces became separated, which in this type of model was not possible because both tissues were supposed to remain in contact. 5.6.13.3 Fibre direction calculation As it was explained in section 5.6.12, a correct modelling of the breast tissue could need allowing for a fibred model. However, as stated in that section, no information was found in the literature about the preferential direction of the fibres in the adipose tissue in general, and in particular in the breast. Here, a method was proposed to estimate the preferential directions of fibres in the adipose tissue based on the works of some authors who observed how living tissues are remodelled to adapt its microstructure in the most efficient way from a mechanical point of view [46, 134, 135, 136]. In particular, tissues reinforced with collagen fibers are thought to have the fibres oriented in the direction in which the tissue is withstanding the largest tensile stress. In this manner, the mechanical properties of collagen fibres are exploited in the most efficient way.