THE EFFECT OF MECHANICAL INTRINSIC FACTORS INDUCED BY MORPHOGENESIS ON BRAIN MECHANICS By Atacan Yücesoy A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering—Doctor of Philosophy 2024 ABSTRACT The brain soft tissue is subject to large strains due to the shape changes that occur during growth. The biological growth alters the state of stress and leads to residual fields existing in the equilibrium state of the tissue after morphogenesis. This resulting stress field typically involves compressive and tensile stresses that vary through the material in a complex fashion. Hence, the mechanical response of residually stressed tissues to finite deformations differs from that of stress-free tissues. Furthermore, considering the role of mechanical properties of the tissues on the regulation of the essential behavior of the cellular structure, the residual fields have a potential role in mechanotransduction at the tissue and cellular scale. The residual fields also should be included when seeking to model the micromechanical mechanisms that give rise to brain injury. While the physical mechanisms of acute and secondary injuries led by extreme events (e.g., blunt impact, blast waves, cavitation, etc.) still remain unclear, it does seem clear that residual stresses could have a significant effect. For example, a preexisting tensile residual stress could accelerate the formation of microfissures during an episode of physical trauma, whereas a preexisting compressive residual stress field could provide some benefit in delaying fissure formation. The research work presented in this dissertation concerns the study of morphogenesis-induced residual stress fields in hyperelastic materials and the potential effect of these residual stress fields on the material response. This research is generally based on the non-linear theory of elasticity. To address the effect of the residual stress field on the material response, the solution of a finite deformation boundary value problem for a residually stressed elastic spherical shell subject to pressure inflation is first provided. To this end, the general constitutive equation for an isotropic Mooney-Rivlin type of hyperelastic material with a background residual stress field is derived. The problem is expressed as a compact integral expression including the base response of the material and the response arising from the presence of the residual stress field. An asymptotic analysis is conducted to examine the dependence of residual-stress integrals on a dimensionless measure of radial strain. The results are compared with the base response of the Mooney-Rivlin type material to pressure inflation and the potential effect of the residual stress field on the material response is discussed. The numerical analysis shows that the residual stress fields have the potential to alter the qualitative behavior of the pressure-inflation response of the material. The work then proceeds with the examination of residual stress due to differential growth in adjoining tissue in incompressible isotropic hyperelastic single and bilayer spherical shells. The kinematics of differential volumetric growth utilizing the incompressible hyperelastic framework are presented for each geometry considered, and the growth-induced residual stress fields are computed for five different growth conditions. Then, the sensitivity of the resultant stress field to the differential growth in adjoining layers is examined for the combination of the five growth conditions. To address the residual stress fields generated by the morphogenesis including symmetry-breaking bifurcation and beyond, the study later continues by building an elementary computational model with idealized geometry, boundary conditions, and parameters. This 2D plane strain computational model provides the residual stress/strain fields emerging in a formation resembling sulcus-gyrus structure in a gyrified brain. Following the differential growth hypothesis, the residual stress and strain fields are computed for the domain where the cortex undergoes only tangential (in-plane) growth while the subcortex does not experience any growth. A detailed stress and strain analysis of the resultant sulcus-gyrus formation is performed to understand morphogenesis-induced residual fields specifically for the sulcal floor and gyral crown. Due to the specific attention to physical injuries leading to the neuropathies such as Chronic traumatic encephalopathy (CTE), the analysis is extended to encompass the response of non- residually stressed sulci subjected to intrasulcal deformations. A 2D plane strain computational model of a single sulcus is built to examine the deformations associated with the expansion of a cavitation bubble in the intrasulcal region. Based on the experimental data, the quasi-static and transient pressure loading conditions are implemented to the gray matter-cerebrospinal fluid (CSF) boundary, and the response of the sulcus is investigated in detail. The findings demonstrate that cavitation result in sulcal expansion and the formation of localized high strain and strain rates at the depth of sulci. The strain and strain rate localization regions resemble the tauopathy / neurofibrillary tangles patterns seen in early CTE. Copyright by ATACAN YÜCESOY 2024 ACKNOWLEDGEMENTS I would like to express my deepest gratitude and appreciation to all those who have contributed to the completion of this doctoral dissertation. Without their support, guidance, and encouragement, this research would not have been possible. I wish to express my profound gratitude to my supervisor, Dr. Thomas J. Pence. His profound expertise, unwavering support, and insightful feedback have played a pivotal role in shaping not only this dissertation but also my overall academic development. I am profoundly indebted to him for his invaluable mentorship and guidance. Additionally, my heartfelt thanks go to my co-advisor, Dr. Ricardo Mejia-Alvarez, whose valuable insights, supports, always constructive critiques, and thoughtful recommendations have significantly enhanced the quality of my work. I would also like to extend my sincere appreciation to the remaining members of my thesis committee: Dr. Adam M. Willis and Dr. Michele J. Grimm. Their devoted time and unwavering commitment to reviewing my research and providing invaluable input have been immensely appreciated. Lastly, I extend my heartfelt appreciation to the faculty and staff of the Department of Mechanical Engineering for their unwavering support and the resources they have provided, which have played a crucial role in the successful completion of this dissertation. I want to extend my heartfelt appreciation to my dear colleagues and friends: Suhas Vidhate, PhD, Joseph Kerwin, PhD, Kian Kalan, PhD, Faezeh Masoomi, Basil Abdelmegied, PhD, Bianca Davila-Montero, and Tushar Kailu. They have been steadfast companions on this demanding journey of mine. Their unwavering encouragement, enlightening conversations, and unwavering moral support have been the guiding stars that kept my motivation burning bright during the toughest of times. I am genuinely thankful for their friendship and the wonderful camaraderie we share. I wish to express my heartfelt gratitude to my family, whose unwavering love, understanding, and unwavering encouragement have been a beacon of light throughout this entire journey. My cherished wife, Sevil Yücesoy, and my precious daughter, Derin Yücesoy, have been steadfast companions, offering their boundless and unconditional support. I consider myself truly blessed to v have them by my side. Lastly, I want to extend my appreciation to my parents Reyhan and Altan Yücesoy, and other family members Aydan, Cumhur and Ece Sultan Özer, Kemal and Kader Kaya, Ulaş, Zeynep and Çınar Ergen for their unconditional support, and love. In conclusion, I am humbled and honored to have had the opportunity to undertake this research, and I recognize the collaborative efforts of all those mentioned above. Their contributions have played an integral role in shaping this dissertation and my academic journey. I am deeply grateful for their support and guidance, and I remain indebted to them for their invaluable contributions. vi TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 GENERAL TREATMENT OF RESIDUAL STRESS IN HYPERELASTIC SPHERICAL SHELLS . . . . . . . . . . . . . . . . . 25 CHAPTER 3 RESIDUAL STRESS FIELDS IN THE UNIFORM GROWING OF A HYPERELASTIC SPHERICAL SHELL . . . . . . . . . . . . . . . . . 52 CHAPTER 4 MICROMECHANIC INVESTIGATION OF THE GYRI-SULCI . . . . 74 CHAPTER 5 THE GROWTH-INDUCED RESIDUAL STRESS FIELDS RESULTING FROM SULCUS-GYRUS FORMATION IN THE BRAIN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 CHAPTER 6 CONCLUDING REMARKS, SPECULATIONS, BROADER CONNECTIONS AND FUTURE DIRECTIONS . . . . . . . . . . . . . 129 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 vii CHAPTER 1 INTRODUCTION 1.1 Introduction and motivation Cortical growth and gyrification lead to a highly convoluted outer surface, a salient hallmark of brain morphology that exhibits a strong correlation with higher-order cognitive capacities and functional intricacy. The cortical folding increases the neocortex volume while minimizing the spatial occupancy of the neuronal network [77]. This structural configuration causes an energet- ically efficient and spatially frugal configuration for the brain. It may, thus, be conceptualized as an optimization process wherein the functionality is maximized, whereas the metabolic cost is minimized within the limited volume. Cortical growth and folding lead to two characteristic structures on the outer surface of the brain: gyri (elevated ridges) and sulci (depressed grooves). Along with the morphological changes, the mechanical state of the tissue also changes due to the slow and large deformations during the morphogenesis. This emerging mechanical state becomes the natural state of the tissue by persisting after the morphogenesis and thereby determining the homeostasis conditions at the tissue and cellular level. In the absence of the surface tractions, the stress field that exists in the equilibrium is known as residual stress [63]. Residual stress fields exhibit both inhomogeneous and anisotropic characteristics. It encompasses compression and tension fields simultaneously due to the zero- traction condition. Moreover, it leads to the presence of localized internal stress fields within the material [75, 74]. Thus, it becomes evident that the mechanical behavior of bodies bearing residual stress is notably more complex compared to the stress-free bodies. Consequently, a profound understanding of the residual mechanical field becomes imperative in order to comprehend the mechanical response of a residually stressed body when subjected to finite deformation and the mechanical homeostatic conditions of tissues. The residual stress often exhibits a complex spatial distribution. Hence, the traditional invasive techniques (e.g., cutting experiments) may not be sufficient to identify the natural stress-free configuration of the tissue completely. Previous experimental investigations have focused on 1 characterizing residual stress fields within various tissues, primarily tubular structures such as arteries [32], aorta [43], and trachea [59]. However, there is a noticeable scarcity of experimental studies aiming at the residual stress field in the cortex. Xu et al. conducted a series of dissection experiments on coronal slices of both the adult mouse brain [143] and developing ferret brain [144] to map the residual stress fields. These experiments involved radial cuts to the tissue with various depths, revealing that the gray matter experiences compression while the white matter undergoes tension [143]. Still, it is important to note that these findings are confined to specific cutting directions and regions. In contrast to tubular tissues, where the opening angle serves as a reliable measure of residual stress, the resulting deformation following cuts in brain slices does not offer a dependable metric for determining the magnitude of residual stress due to the brain’s intricate geometry. The findings also suggested that the observed residual stress field is likely a result of differential growth. The differential growth hypothesis has been tested in experiments involving swelling layered gel samples, which have led to the development of brain-like morphologies [125, 126, 58]. From a mechanical point of view, there have been several hypotheses in the literature explaining cortical growth and folding, which have been recently reviewed in [10, 122]. These hypotheses encompass the following key facets: a) constrained cortical expansion, b) axonal tension, and c) differential tangential expansion. The details of each hypothesis will be discussed in detail later. The constrained cortical expansion hypothesis suggests that the external constraints generated by the surrounding structures (e.g., skull) lead to compressive stress on the outer cortex surface, which shapes the cortex. Yet, the experimental evidence showed that cortical folding can take place in the absence of the outer constraint [8]. The axonal tensional hypothesis suggests that tension generated by growing axons plays a significant role in shaping the morphology of the cortex, particularly in the formation of gyri and sulci. However, relatively recent experimental studies showed that the axonal tension is not sufficient to form the gyri and sulci. Furthermore, the tensional stress mostly exists in the subcortical white matter tract instead of the core of the developing gyri [144]. The differential hypothesis basically asserts that the cortex grows more than the subcortex in the 2 tangential direction. The residual stress field observed in the cutting experiment is consistent with the differential cortical growth hypothesis [144]. Considering the mechanical aspect of the growth mechanism, the effect of mechanical parameters on growth-induced folding needs to be addressed. Noteworthy among these are the potential influence of stiffness disparities between gray and white matter, the initial thickness of gray matter, the fetal brain’s configuration at the onset of folding, and the growth ratios across distinct layers. These quantities have been extensively investigated in [21] and [140]. In terms of the mechanics of gyrification, the aforementioned works in the context of biological growth are mechanically similar to surface instabilities in other elastic systems. This includes multilayer films subject either to a strain mismatch due to various loading stimuli (e.g., growth, swelling, atrophy, thermal expansion, initial stress, other chemical and physical processes of film deposition, etc.,) that can lead to compression in the outer layer [13, 4]. Irrespective of the cause of the compression, the generic behavior of such a system is to respond to a mismatched strain through the formation of surface wrinkles. As the mismatch strain increases, the progression of surface instability in the system beyond the wrinkled state depends on the stiffness contrast and interfacial energy [137]. For a small stiffness ratio, the system develops more localized deformations such as creases and foldings [81]. The intricate convolutions observed in the brain’s development are thought to arise, at least in part, as a consequence of these instabilities [51, 58]. For larger stiffness ratios, the system exhibits more intricate phenomena, including period doubling and ridge-type instabilities [23]. On the other hand, there is a substantial body of computational research dedicated to com- prehending the mechanisms driving cortical growth, owing to its significant clinical implications. To this end, both elementary and advanced computational models of cortical folding, known as gyrification, have been employed to investigate various aspects, encompassing (i) the replication of observed cortical folding patterns in actual gyrated brains [125, 126, 133], (ii) the detection of anomalies associated with cortical development ([19, 135]), (iii) the fundamental processes gov- erning cortical growth and folding [152, 71, 90, 21, 139, 125, 126], and (iv) the evaluation of the 3 influence of physical and material parameters, alongside their regional variations [140, 33, 21, 20]. A comprehensive overview of the common methodologies employed in numerical investigations of cortical folding can be found in a recent review by Darayi et al. [34]. Despite the advancements in these models in terms of complexity and predictive capabilities, their ability to discern residual stress patterns during brain development remains somewhat limited. Consequently, addressing residual stress fields is imperative due to their potential functional and mechanical implications. Bayly et al. [9] developed a numerical model to elucidate stress patterns in the developing brain, although their findings primarily pertained to sinusoidal wrinkling and isolated gyrus formations. Furthermore, [125] identified the presence of residual stress fields in cusped foldings but did not extensively delve into their distinctive characteristics. In this study, a numerical analysis is conducted to examine the following aspects: (i) the potential effect of residual stress field on the mechanical response of the isotropic hyperelastic soft material to pressure-inflation deformation, (ii) the attributes of the residual stress field induced by volumetric growth in idealized geometries (e.g., single layer and bilayer spherical shells), (iii) the strain/stress fields generated by the cavitation-induced damages occurred in the single sulci, (iv) the residual stress fields emerged in the soft material due to the inhomogeneous growth including the symmetry-breaking regimes. In the following part of this chapter, an extensive review of the literature concerning the mechanical aspects of cortical folding and finite volumetric growth is presented. It is followed by the review of a comprehensive examination of computational models that explore cortical growth and gyrification from various aspects and the experimental works that investigate the behavior of residual stress fields within the cortex. In the rest of this chapter, we introduce a somewhat broader theoretical framework to incorporate the kinematics of volumetric growth. In Chapter II, we provide a solution for the boundary value problem for a residually stressed elastic hollow spherical shell subject to pressure inflation by employing a conventional hyperelastic 4 model. Specifically, we consider a well-established Mooney-Rivlin-type material model, whose response to pressure inflation is characterized by the invariants and relative shell thickness. In the analysis, the reference configuration is considered to be inherently residually stressed. To this aim, four simple different residual stress profiles with various levels of strength are utilized. The strain energy function is redefined to incorporate the invariants resulting from the residual stress field considered in addition to the finite deformations. After presenting the boundary value problem and a comprehensive integration procedure in detail, (i) an asymptotic analysis is conducted to examine the sensitivity of the residual stress integrals on the radial strain, (ii) the potential effects of the residual stress field on the baseline pressure-inflation response of the material is examined by numerical analysis. In Chapter III, the growth-induced residual stress fields in an idealized single and bilayer incompressible hyperelastic spherical shell are examined using the nonlinear theory of elasticity. In this regard, the radial and circumferential residual stress fields are computed for five different growth conditions: area, surface, isotropic, and the combination of area and surface. A numerical analysis is conducted to understand the effect of the differential growth conditions in the bilayer spherical shell. Lastly, the numerical results are discussed and compared to previous experimental findings. In Chapter IV, to address cavitation-induced damages specifically occurring in the intrasulcal regions and its possible link to neuropathies/tauopathies such as CTE, an elementary 2D plane strain computational model is built using small deformation theory. In the computational model, the response of the sulci under quasi-static and time-dependent pressure boundary conditions is examined by simulating the deformation generated by the expansion of the cavitation-induced bubble in the intrasulcal region. Numerical results predict highly localized strain and strain rate concentrations especially at the deepest point of the sulcus due to the expansion, which can lead to a neuropathology like CTE. In Chapter V, the study later continues with the analysis of the growth-induced residual stress fields comprising the effect of the cortical folding process. The primary objective of this chapter 5 is to investigate the resulting morphology and mechanical state of the tissue, particularly within folding pattern that resemble the sulcus-gyrus formation. To this purpose, an elementary 2D plane-strain finite element model is built to compute the displacement field leading to the lowest total storage energy of the system among all kinematically admissible states. The computational domain is modeled as initially flat bilayer rectangle, and adopts the differential tangential growth hypothesis. The differential growth rate between the cortex and subcortex gives rise to a strain mismatch, serving as the primary parameter for distinguishing various grown states presented throughout this chapter. The chapter then proceeds with presenting the morphological changes in evolving folding patterns and the residual stress/strain field for critical mismatched strain states. For the mature sulcus-gyrus formation, the analysis is proceed with the computation of the eigenvalues and eigenvectors of the stress field, enabling an investigation into the principal stress fields and their associated directions. Lastly, the effect of skull constraint on the folding pattern and its mechanical state is examined. 1.2 Literature review 1.2.1 On the mechanics of cortical growth and gyrification Early attempts to explain the underlying mechanics of cortical folding posited that local dispro- portions in the rate of size change, resulting from the molecular and cellular transformations during cortex development, were responsible for cortical folding. Smith et al. [120] proposed that cortical folding arises from the disproportion between the volume of the cortex and the surface area of the neocortex, rather than being solely attributed to local changes. On the other hand, subsequent analytical models primarily focused on extrinsic factors that could shape the cortex. One such model postulated that constraints imposed by the surrounding brain structures, such as the skull, induce buckling in the expanding cortex’s outer surface [80]. In other words, the relatively rigid restrictive structures such as the skull could generate compressive stress on the growing cortex, leading to buckling. In 1950, Barron [8] experimentally examined the effect of the skull on a growing cortex using surgically manipulated fetal sheep brains. A part of the left hemisphere of the fetal brain in utero was removed to reduce the pressure generated by the skull and to provide more 6 intracranial space for the intact right hemisphere. Contrary to earlier hypotheses, the experimental findings demonstrated that (i) the external constraint is not necessary to initiate the cortical growth and (ii) a highly convoluted outer surface still developed despite the surgical intervention, with only minor morphological differences observed. It is fair to assume that external constraints play a role in shaping the final brain morphology, such as flattening gyral crowns, even if they do not directly initiate cortical folding. Subsequently, Richman et al. [106] introduced a mechanical model based on the differential elasticity between the inner and outer stratum. In this model, the cerebral cortex was assumed as a relatively thin and stiff layer (gray matter) growing over a relatively soft infinitely-deep elastic core (white matter). The model proposes that the outer layer growing more than the core leads to a sufficient compressive stress field within the outer layer to initiate the cortical buckling. This model is capable in predicting the wavelength of the foldings using the elastic properties of the tissues and thickness of the layers. However, this model had limitations in two aspects: (i) the assumptions for the relative elastic properties of the gray and white matter ranging between 1/100 to 1/10. However, these properties does not reflect the actual material properties. (ii) Due to the elastic properties, the model estimated a sinusoidal wrinkling pattern that does not reflect the actual brain folding pattern. Currently, there are two prominent hypotheses pointing out the driven mechanism associated with cortical folding during the evolution of the brain as illustrated in figure 1.1. One is axonal tension and the other is the differential growth hypothesis. The axonal tension hypothesis states that the tension in the axons pulls the interconnected cortex regions, thereby reducing the spatial separation between opposite banks of gyri [131]. However, recent experimental studies showed that the residual stress field observed in an evolving ferret brain is consistent with the differential growth hypothesis, and the axonal tension is not sufficient to generate a highly convoluted cortex as seen in the mammalian brains [144]. On the other hand, multiple studies have indicated that axons can influence cortical folding patterns beyond the generation of tension [122]. An experimental study has reported deviations in the folding morphology of the occipital cortex in monkeys whose eyes were surgically removed during embryogenesis, resulting in patterns different from the expected 7 pattern [101, 36]. Furthermore, the reduction in size and number of axons and nuclei associated with sight has been observed. A relatively recent experimental study indicated that the dissection of eyes reduces the size of cortical gyrus in the primary visual cortex [105]. Hence, while axonal tension may not serve as the primary mechanism, it undeniably plays a significant role in shaping the cortex. Figure 1.1 Illustration of cortical folding mechanisms: axonal tension (left), differential tangential growth (right). The kinematic description of volumetric finite growth was first formulated by Skalak [117], and his coworkers [118] in the context of continuum mechanics. Following these early attempts, Rodriguez and Hoger [108] presented the general formulation for the finite kinematics of the volumetric growth, further expanded upon by Chen and Hoger [30]. The growth of a biological tissue is described as two distinct phases: i) the addition of tissue material due to the inward mass flux (e.g., cell migration), or the mass sources (e.g., cell growth, division, or enlargement) and ii) subsequent to growth, an elastic deformation is required to restore tissue compatibility. It is important to note that the growth field need not be continuous, which can result in discontinuities within the body. Therefore, an elastic deformation is necessary to make the body compatible. Furthermore, these elastic deformations can induce a residual stress field within the body, which is initially assumed to be locally stress-free prior to any growth. These two subsequent processes are mathematically modeled as the multiplicative decomposition of two tensor valued functions. 8 From a mechanical perspective, the mechanics governing cortical growth and folding exhibit similarities to surface instabilities seen in other highly elastic systems. Surface instabilities can emerge in either homogeneous [83] or multilayered [136] highly elastic soft materials when sub- ject to compressive stress. Various loading stimuli such as growth, swelling, atrophy, thermal expansion/shrinkage, initial stress, and other chemical and physical processes associated with film deposition, can induce compressive strain triggering surface instabilities. The instability conditions for the highly localized buckling deformations in layered viscoelastic and linearly elastic medium date back, which were studied in [12] and [14], respectively. The onset of buckling conditions for a symmetric hyperelastic three-ply construction was examined by Pence and Song [98] in the context of the nonlinear theory of finite elasticity. This study showed the possibility of two types of buckling deformations in this specific case: barreling and flexural buckling. The examination was further extended to address the family of buckled configurations for two-ply laminate associated with the wrinkling load [148]. Regardless of the specific cause of the compressive stress, the typical behavior of such bilayer elastic systems is to respond to the mismatched strain by forming periodic wrinkles on the surface. Assuming a thin and relatively stiff film upon a much thicker and more compliant substrate, the initial wavelength of the wrinkle can be estimated utilizing the formula 𝜆 = 2𝜋𝑡film(𝜇film/3𝜇sub)1/3, (1.2.1) where 𝜇film and 𝜇sub represent the respective stiffness modulus of the film and substrate, while 𝑡film denotes the thickness of the film [29]. The elastic system releases the compression of the film via periodic wrinkles to lower the total energy of the system, because bending is energetically less costly than compression. The wavelength of the initial wrinkles can be predicted based on linear stability analysis [23]. With further increases in the mismatched strain, the wrinkles become increasingly unstable, giving rise to more complex surface patterns such as folds [100], ridges [23], and period-doubling [15], depending on the elastic properties of the multilayered system. Alternatively, an initially flat surface in the elastic system may transform directly into a crease [136], resulting in localized and discontinuous surface transformation without the formation of wrinkles [72]. It should be noted that sulcification, as suggested by Hohlfeld and Mahadevan [65], 9 represents a scale-free nonlinear subcritical instability distinct from wrinkling and buckling. Wang et al. [136] introduced a diagram to distinguish various growth-induced surface instability patterns based on the mismatch strain, stiffness ratio, and normalized adhesion energy between the layers for a bilayer structure. Soft materials are prone to surface instabilities triggered by buckling due to the low elastic moduli and being highly responsive to external stimuli. Therefore there has been an increasing interest in surface instabilities because of their applicability to soft tissues [81]. In the context of cortical folding, the gyrus and sulcus formations can be seen as surface instability arising from the constrained differential growth . The difference in growth rates of cortex and subcortex leads to compressive mismatch strain in the outer layer that leads to an instability when exceeding the certain threshold. Swelling gel experiments, which replicate the volumetric expansion of the cortex through constrained differential growth, demonstrate how surface instability patterns evolve as the mismatch strain increases [58, 125], leading to a convoluted surface that resembles a human brain [126]. In [58], it is demonstrated that cortical folding likely results from a wrinkling instability followed by a folding transformation, rather than a direct transition from a flat to a creased state. Ben Amar and Goriely [6] incorporated the finite growth into the incremental deformation theory and showed how an anisotropic growth field causes surface instabilities in an incompressible hyperelastic spherical shell by means of affecting the geometry and the internal stress field of the shell. Assuming the gray and white matter is perfectly bonded, the diagram introduced by [136] predicts a crease type instability for a low stiffness contrast, i.e., 𝐸gray 𝐸white ≤ 1.3, while wrinkling and folding arises for the stiffness ratios 𝐸gray 𝐸white > 1.3 for the primary bifurcation. As the mismatch strain increases, more advanced modes of instabilities such as period-doubling and period-tripling can emerge depending on the stiffness ratio in the post-bifurcation regime [18]. 1.2.2 Computational Studies on Cortical Growth and Gyrification A variety of motivations have turned into elementary and advanced computational models aimed at simulating the mechanics of cortical growth and folding to gain insight into (i) the formation of the cortical folding patterns as seen in the actual gyrified brain [125, 126, 133], (ii) developmental 10 abnormalities linked to cortical growth and gyrification [19, 135], (iii) fundamental mechanism of the cortical growth and folding [21, 71, 90, 125, 126, 129, 139, 152], and (iv) the potential effects of physical and material parameters (e.g., the initial thickness of the cortex, stiffness ratio, and growth rate difference of the layers, and initial shape of the cortex) on the folding patterns [21, 20, 33, 140]. The computational studies in the field of cortical growth and folding has been reviewed in [35]. Continuous and discrete mathematical models are two approaches to simulating biological growth and remodeling processes. Discrete models are primarily employed to replicate biological growth at discrete levels, particularly at the cellular and subcellular scales, (e.g., tumor growth). Conversely, continuum modeling operates at the tissue scale, treating the tissue as a continuous medium during the growth process. The choice between continuous and discrete approaches represents an ongoing research aspect within the field of developmental biology [11]. In the context of brain soft tissue growth and folding, the continuum approach is mostly used. The theory of finite volumetric growth, introduced by Rodriguez and Hoger [109], is the current approach used for describing the kinematics of biological growth and remodeling. This kinematic model combines the effect of growth and elastic deformation, which is represented by a multiplicative decomposition of the deformation gradient tensor. The growth of the cortex is defined by a model aligning with the differential growth hypothesis, which is more commonly used compared to the other mechanic models presented in Section (1.2.1). Experimental study shows the residual stress fields observed in the ferret brain consistent with the differential growth hypothesis [144]. This has been later supported by the computational models [125]. The differential growth basically proposes that the gray matter grows more than the white matter in the tangential direction. This requires an individual tensor-valued function for the growth of each layer in the model. The components of the growth tensor could be a function depending on such as the spatial location [125], the growth of the adjacent layer [19], or an independent function such as logistic-growth function [90, 129]. The mechanics of cortical growth and gyrification can be also modeled by coupling with differential equations representing, at least partly, the effect of underlying biological changes. Verner et al. [133] introduced a cortical growth and folding model coupled with the 11 advection-diffusion-reaction equation to include migration and tangential intercalation of neuron cells during the development. On the other hand, an alternative approach considers the growth rate as a singular scalar value, a method employed in several studies [19, 104, 135]. This approach implies isotropic growth for the cortex. In cases where the growth rate is constant, it must exceed one to signify volumetric positive growth. A growth rate that equals unity indicates no growth, while values below the unity represent atrophy and resorption. Furthermore, there are alternative approaches to replicate the mechanics of the volumetric expansion of the tissue during the growth and morphogenesis such as thermal expansion [71, 135], and osmotic expansion [46]. The volumetric expansion is controlled by the temperature of the layer and osmotic pressure respectively. These factors act as the growth functions instead. One of the advantages is that such a treatment allows one to make use of commercial finite element codes (e.g., ABAQUS, Comsol) because of the availability of built-in algorithms. To the best of my knowledge, the axonal tensional hypothesis is employed in only one work [46]. The axonal tension was modeled by applying tension forces along the white matter fiber orientation measured by diffusion tensor imaging (DTI). Because of the lack of experimental data related to the tensile force in the axons, the magnitude of the tensional force was determined in the same order as the cerebrospinal fluid (CSF) pressure applied to the cortex. The results showed that the differential tangential growth leads to a sulcal pattern resembling the actual human brain rather than the axonal tension. The shape and geometrical attributes of the cortex prior to any growth are effective on the resultant folding patterns and morphology, which is also a research question. In computational studies, the cortex has been modeled using (i) various simplified geometries representing a certain subsection of the brain or whole and (ii) the geometries derived from Magnetic Resonance Imaging (MRI) images. The cortex has been extensively represented using various simplified geometries either in two- or three- dimensions such as a semi-sphere [133], spherical shells [125], ellipsoid [9, 135, 140], rectangular blocks [19, 66, 71, 135, 139], circular disk [125, 129, 144], coronal slice [143], semi-ellipse [71], and hemispherical [133]. The rectangular domain with plane strain 12 assumption is the most used model in the literature [19, 20, 33, 70, 104, 125, 139]. In addition to the idealized geometries, MRI-based models of the brain, which have a complex shape even before folding, have been also used to derive more realistic reference configurations. Of the reviewed studies, there are two FE models where MRI images of a fetal brain at 22 weeks of gestation are employed to obtain a more realistic reference configuration prior to any folding [90, 126]. In all reviewed studies, the brain is treated as a bilayer structure that consists of an outer layer with various thicknesses on a relatively thick substrate. In the context of the brain anatomy, the outer layer and substrate correspond to the gray and white matter respectively. In the computational models, it is commonly assumed that the gray matter thickness is uniform across the entire domain [19, 71, 135]. It should also be noted that the initial thickness of gray matter is one of the parameters that have been extensively examined in the models. On the other hand, there have been computational models where the thickness of gray matter exhibits heterogeneity along the white matter [33]. Unlike the variety in the geometry of reference configuration, it is fair to say that there is a uniformity in the selection of constitutive form to model the brain soft tissue. Considering the time scale of a typical biological growth process, it is fair to assume that the time scale of growth is acting longer than the other time scales (e.g., stress relaxation, viscous and elastic time scales) [50]. Hence, cortical growth and folding is generally considered a quasi-static process. That leads to a tendency to a time-independent constitutive form. Along with the linear elastic model, the hyperelastic constitutive models have also been used for the brain soft tissue. More specifically, the well-known neo-Hookean type strain energy function has been heavily used in the computational model [19, 126, 104, 103, 138]. Such a model could be sufficient to estimate how the gyrification happens without requiring a high computational resource. The quantitative values of the material properties are not discussed in this review purposely. Barron’s experimental findings [8] established that even if cranial confinement is not required for the initiation of cortical folding, the cranial constraints can exert (in part at least) influence, particularly on the ultimate cortex morphology. The presence of relatively stiffer surrounding tissues and skulls may contribute to the development of flattened gyral crowns [122]. Certain 13 computational models have included the influence of the skull on cortical folding patterns. In the literature, this volumetric constraint has been simulated through two primary approaches: (a) the inclusion of a rigid shell to represent the skull or (b) the application of displacement constraints on nodes. Tallinen et al. [125] noted that cranial confinement leads to the formation of flattened gyral crowns. Nie et al. [90] suggested that the skull has an important regulatory role in cortical folding even if the skull constraint is not necessary to trigger the folding mechanism. The results show that skull constraint results in more intricate folding patterns, deeper sulci, and earlier development of primary and secondary sulci compared to the non-skull-constrained models. Cerebrospinal fluid (CSF) could be another source that can generate a constraint on a developing brain. The CSF applies a pressure ranging from 1.5 mmHg to 6 mmHg (corresponding to 200 Pa to 800 Pa) in the case of a full-term infant brain and from 10 mmHg to 15 mmHg (equivalent to 1300 Pa to 2000 Pa) in the adult brain [38]. A recent computational investigation conducted by Jafarabadi et al. [70] has explored the effect of cerebrospinal fluid (CSF) on the buckling behavior of a bilayer system. The results show that the presence of CSF has an impact on the stability of the bilayer system, thereby the buckling behavior. Still, the attention on the role of CSF is scarce. The layers of meninges (the pia mater, dura mater, and arachnoid mater) have the potential to introduce constraints during cortical growth and folding processes. However, any studies that have incorporated the consideration of meninges layers within computational models do not exist in the literature. In addition to these considerations, reviewed computational models typically employ two com- mon boundary conditions: (a) perfect bonding between the cortex and subcortex to ensure interface compatibility and (b) prevention of self-contact between the outer layers of gray matter. Certain studies also introduce a cortical heterogeneity that acts as a small perturbation. These perturbations have the potential to trigger the formation of cortical folding. Such perturbations may be in the form of a localized variation in geometrical and material attributes such as thickness [33], stiffness, cortical growth rate [20], or the application of external mechanical forces [125]. It is noteworthy that such perturbations (e.g., inhomogeneities in thickness, stiffness, and growth rates) influence the localization of specific sulcus-gyrus formations, and the growth-induced instability patterns, 14 which also leads to more complex formations [20, 33]. Previous studies revealed that gyri tend to develop more frequently in regions where the growth rate and stiffness are high, while sulci tend to emerge at locations characterized by lower growth or stiffness [9, 20, 104]. Heterogeneous thickness distribution can introduce mechanical perturbations that affect the initiation of folds. In other words, folding initiates at locations where the top layer is relatively thinner, once sufficient compressive stress accumulates in the outer layer [33]. The growth-induced cortical folding pattern is significantly influenced by several factors. These factors include a) Cortex Thickness: A thicker cortex results in a reduced number of gyri and an increased presence of isolated sulci. In essence, as the cortex thickness increases, the wavelength of instability at the free surface also increases [71, 140, 19, 133]. b) Stiffness Contrast: A higher stiffness contrast between the cortex and subcortex leads to an increased wavelength, with wrinkles becoming the dominant type of surface instability [71]. c) Growth Rate: Higher growth rates are associated with the formation of deeper sulci and smaller gyri [90, 135], and the emergence of gyri is more likely in regions characterized by higher growth rates [20, 152]. Numerous validation approaches have been implemented to assess the degree of proximity between computed results and the actual brain surface. These validation methods predominantly rely on qualitative assessments [19, 135], the analysis of residual stress fields within the formations [125], morphometric techniques [126], as well as the utilization of established metrics such as the gyrification index (GI) [19], and novel metrics [140] to quantitatively characterize features of the resulting folding pattern. One commonly utilized quantitative measure for assessing folding patterns is the sulcification wavelength, defined as the distance between two adjacent gyri. Additionally, other quantitative metrics include dimensionless mean curvature, the surface-based 3D gyrification index, and sulcal depth [140]. Wang et al. [140] also introduced a novel metric designed to quantify the anisotropy of folding orientation, facilitating an exploration of the relationship between fold orientation and brain shape. Tallinen et al. [125] conducted a swelling gel experiment aimed at capturing qualitative characteristics of the folding pattern, providing support for the underlying physical mechanisms governing cortical folding. Their work also posited that the intricate surface 15 morphology of the human brain is an inevitable consequence arising from the constraints imposed by differential growth [126]. 1.2.3 Experimental Studies on Residual Stress Fields in Brain The internal histological changes associated with the external features of developing gyrus have been examined to understand how a gyrus and sulcus form during morphogenesis. Thus, this approach allows us to make a connection between the deformation aroused by the growth and the internal changes in the tissue. In an experimental work conducted by Smart and McSherry [119] performed on postnatal ferret brain, a series of changes has been observed: a) the gyrus is formed as the result of the radial movement of the cortical plate relative to the sulcal floor, b) the radial tissue lines tend to be curved instead of increasing the length during the sulcus formation. It also showed that the movement of cell nuclei is tangential in the sulcal area rather than the radial [62]. It is worth noting that the ferret has a lightly folded brain compared to humans. In that regard, further analysis is still require to understand the advanced foldings seen in the human brain. Additionally, there are geometrical models proposed that the initial geometry of the cortex could affect the pattern of folding. The general direction of folding follow the minimum curvature of the brain considering a uniform brain growth and minimal lateral movement between the outer layer and inner core [128]. This model is also consistent with the growth-driven model proposed by Richman. The presence of residual stress in soft biological tissues can be demonstrated through an invasive experimental technique. This method consists of the incision of tissue in various orientations, causing the tissue to change shape as it relieves elastic stress. Subsequently, one can compute the residual stress within the material by tracking the alterations in shape and obtain a constitutive law using this stress-free configuration. This approach has been employed to detect the residual stress fields, mostly the tubular structures such as the arteries [134], ventricle [92], ureter [57], and trachea [59]. For more details about this invasive method and its applications, I refer to [44] and the related chapters of [45]. The behavior of the incision allows for predictions regarding the characteristic (tensile or compressive) of the residual stress in the perpendicular direction to the cut, depending on whether the incision remains closed or not. The size of the opening is utilized as 16 a quantitative measure to evaluate the magnitude of residual stress. Yet, it should be noted that the size of the opening might not be a precise metric for quantitative analysis due to the persistence of a substantial residual stress component even after a singular incision [130]. Moreover, the opening may be influenced by the specific treatments applied to the tissue prior to the cut. Vossoughi et al. [134] showed an increase in the opening angle after making a circumferential cut on the arterial segments that already have been cut in the radial direction. In another experimental study, Greenwald et al. [56] observed variations in the opening angle of incised bovine carotid arteries depending on prior removed the amount of wall material. Returning our attention to the dissection experiments conducted on the brain tissue, it is noteworthy to mention that there have been only two experimental studies in the literature. These studies have been conducted on the adult mouse [143] and developing ferret brain [144]. One of the primary distinctions between these studies pertains to the morphological attributes of the brain specimens in question. Specifically, the mouse brain exhibits a smooth outer surface, whereas ferrets possess a lightly folded cortex in comparison to other mammals, such as primates. With this regard, these works give us an opportunity to assess the residual stress fields within the tissue after growth and folding separately. It is also worth noting that the ferret brain finds extensive utility in the field of traumatic brain injury and recovery due to its morphological characteristics, which more resemble those of the human brain [112]. Xu et al. [143] performed a series of opening angle experiments utilizing dissected coronal slices from the adult mouse brain. The primary objective of their investigation was to examine the existence of residual stresses and elucidate their distribution within the tissue. It is worth noting that, for the purposes of this section, our focus shall be confined to the experimental findings, notwithstanding the presence of numerical results within the study. In the course of these experiments, a series of incisions were made solely in the radial direction, traversing both the cortical gray and white matter at varying depths. The outcomes of these incisions revealed distinct behaviors in the cortical gray and white matter. Specifically, the cortical gray matter exhibited a tendency to remain closed, whereas the cortical white matter displayed an opening response. 17 This divergence in behavior indicated the presence of a compressive residual stress field within the gray matter and a tensile residual stress field within the white matter. The experimental results are depicted in Figure 1.2, which also underscores that the growth process itself can engender a residual stress field, independently of gyrification. It also shows that the growth itself can lead to the residual stress field without the gyrification. Furthermore, it is noteworthy that the distribution of residual stress is not uniform along the incision. Qualitatively, a broader opening was observed at the interface between the corpus callosum and the cortical gray matter. It is essential to underscore that this study offers insights solely into the radial direction of residual stress fields, as no incisions were made in the circumferential direction. In the other study, Xu et al. [144] carried out a series of opening angle experiments on the developing ferret brain, employing various cutting directions to discern the presence and charac- teristics of residual stress fields. To investigate the developmental progression of residual stress fields during growth and morphogenesis, these experiments were conducted at multiple stages of ferret development, specifically on postnatal days (P) -6, -11, -15, -18, and in adulthood (9- and 12-month old). The research outcomes were primarily focused on the P6, P18, and adult ferret brains. This investigation yields valuable insights into two key aspects: (1) the differentiation in the accumulation of residual stress fields within sulci and gyri and (2) the alterations in residual stress fields during the process of gyrification. The experimental results reveal that circumferential cuts made within both the sulci and gyri primarily remained closed, with minimal radial tension observed. Conversely, other cuts in these regions did not exhibit any openings. Radial incisions through sulci induced openings in both the cortex and subcortex. A noteworthy distinction was observed, as the outermost surface exclusively exhibited an opening at P6, indicating the devel- opment of compressive residual stress as the cortex folded. In the case of radial cuts within gyri, openings were contingent upon the depth of the incisions. Specifically, deeper cuts, encompassing the sub-plate of the cortex, resulted in openings; while shallow incisions that penetrated only the gray matter remained closed. These findings collectively imply the circumferential residual stress is formed at the pial surface (the interface between gray matter and cerebrospinal fluid) and within 18 Figure 1.2 The result of a cutting experiment conducted on the mouse cortex. The coronal section of the mouse brain is intact (a1,b1). The radial cuts have been made only through gray matter (a2,b2), deeper white matter tract (a3) and thalamus (b3), and after waiting 15 min (b2-b4). The numerical model of corresponding cuts and the contour of circumferential cuts (c1-c4)(Image sourced from [143], reproduced with permission from SNCSC). the deep layers of white matter following gyrification. The experimental observations align with the differential growth hypothesis. 19 Figure 1.3 Snapshots of microdissection experiments conducted on the coronal ferret brain sections at different postnatal days: P6; P18; and adult (9- and 12-month old) (Image sourced from [144], used with permission under copyright license). 20 1.3 Hyperelastic treatment of growth The kinematics of finite volumetric growth takes into account the combined effect of growth and elastic deformation. This theoretical framework was originally proposed in [109], wherein elastic growth is conceptualized as a sequential two-step process. Consider a stress-free unloaded reference configuration denoted by B0 in which X is a position vector of a generic material point prior to any deformation including growth. Growth can lead to a change in shape and size motivated by various biological mechanisms that are regarded as given for the present treatment. The possible biological mechanisms could be cell division, cell migration, or cell enlargement; however, it is beyond the scope of this work. The growth is kinematically described by a tensor valued function Fg mapping the material point in the reference configuration X into a virtual intermediate configuration B𝜉 for the same material point, i.e., Fg : B0 → B𝜉. It should be noted that biological growth is a complex and tissue-specific process. To determine the growth tensor Fg can therefore be challenging using experimental techniques due to the in- compatible displacement field of the intermediate configuration [28]. Mathematically, the simplest form of growth tensor is expressed as Fg = 𝑐I where 𝑐 is a scalar growth multiplier and I is a second-order identity tensor. The growth multiplier 𝑐 must be positive and greater than the unity to represent the growth. The scalar growth multiplier has been extensively used to simulate cortical growth [19, 66, 103, 135]. The growth induces alterations in the size or relative position of the material points within a body. The volume change resulting from growth can also exhibit spatial variation, i.e., 𝑣 = 𝑣(X), which does not need to be continuous across the body. Consequently, this situation introduces a virtual intermediate state, wherein the body’s compatibility is not maintained after growth. The virtual intermediate configuration needs an elastic deformation to restore the compatibility of the body, leading to the development of elastic strain within the body. Kinematically, the growth tensor is accompanied by an elastic deformation tensor Fe that maps the same material point in the virtual intermediate configuration to deformed configuration x, i.e., Fe : B𝜉 → B. Equivalently, the gradient of the map between the reference configuration and the deformed 21 configuration x = 𝜒(X) is the tensor F = 𝜕x 𝜕X , (1.3.1) where F is the total deformation gradients including both the change of shape and volume at all points due to growth and the elastic deformation restoring the compatibility, as depicted in figure 1.4. Hence, the total deformation gradient F is decomposed multiplicatively into F = FeFg. (1.3.2) As the reader is probably aware, a multiplicative decomposition of the form given in (1.3.2) is commonly used in continuum mechanics to describe a variety of physical processes. Similarly, the Jacobian 𝐽 is also multiplicatively decomposed into the elastic volume change 𝐽 𝑒 = det F𝑒 and growth-induced volume change 𝐽𝑔 = det F𝑔. Note that Fg and Fe have positive determinants i.e., 𝐽 𝑒 = det F𝑒 > 0, and 𝐽𝑔 = det F𝑔 > 0. Further, it is often presumed that the Figure 1.4 Standard schematic representation of the kinematics of finite growth as presented in many sources [49, 67]. elastic accommodation is isochoric, which gives the local constraint equation detFe = 1. (1.3.3) Under this additional presumption, the change in volume after the growth and the volume constraint are given by detFg = 𝑣 ⇒ detF = 𝑣. (1.3.4) 22 The energetic framework of hyperelastic materials is represented by a stored energy density 𝑊 as a function of Fe. The strain energy density is the function of the elastic deformation tensor Fe only via the right Cauchy-Green tensor Ce due to the frame invariance of strain energy density i.e., 𝑊 = 𝑊 (Ce, X) where Ce is the right Cauchy-Green tensor associated with the elastic accommodation Fe, namely Ce = FT e Fe. (1.3.5) The elastic deformation tensor develops elastic strains to make the body compatible subsequent to growth. Under some circumstances that will be explained later, the elastic deformation tensor growth may alter the internal stress field of the body. This can give rise to the introduction of residual stress fields induced by growth. The surface tractions t on the boundary 𝜕B of the deformed body B is defined via 𝝈 and outward unit normal n, namely t = 𝝈n on 𝜕B, (1.3.6) where 𝝈 denotes the Cauchy stress. Note that (1.3.6) is also valid for internal surfaces such as tissue interfaces. In the context of morphogenesis, the growth process occurs on a relatively long time scale ranging from days to months, which is larger than any other time scales (e.g., relaxation time scale of soft tissue). Hence, the process can be considered quasi-static. With the absence of the body forces, the equilibrium equation thereby reduces to div𝝈 = 0, (1.3.7) where div is the divergence operation with respect to current location x. The Cauchy stress 𝝈 is given by 𝝈 = 2 𝐽𝑒 Fe 𝜕𝑊 𝜕Ce . FT e (1.3.8) Throughout the presented study, the attention is confined to isotropic materials. The strain energy function 𝑊 depend upon Ce only through its principal scalar invariants, 𝐼1 = tr Ce, 𝐼2 = 1 2 [(tr Ce)2 − tr(Ce)2)], 𝐼3 = det Ce. (1.3.9) 23 The Cauchy stress tensor given in (1.3.8) turns into 𝝈 = 2 𝐽𝑒 (cid:18) 𝜕𝑊 𝜕𝐼 𝑒 1 + (cid:19) 𝐼 𝑒 1 𝜕𝑊 𝜕𝐼 𝑒 2 Be + (cid:19) 2 𝐽𝑒 (cid:18) 𝜕𝑊 𝜕𝐼 𝑒 2 B2 e + 𝜕𝑊 𝜕𝐽𝑒 I, Be = FeF𝑇 e , (1.3.10) where Be is the left Cauchy-Green deformation tensor. Equation (5.1.2) describes how stress in a hyperelastic material will naturally arise from a growth process followed by elastic accommodation when that growth process is either nonuniform or anisotropic or both. Note that when (1.3.3) holds that the third invariant is equal to one. In this case the strain energy function becomes 𝑊 = 𝑊 (𝐼1, 𝐼2) and (5.1.2) turns into 𝝈 = 2 (cid:18) 𝜕𝑊 𝜕𝐼1 + 𝐼1 (cid:19) 𝜕𝑊 𝜕𝐼2 Be − 2 𝜕𝑊 𝜕𝐼2 e − 𝑝I, B2 (1.3.11) where 𝑝 is a hydrostatic pressure associated with the incompressibility constraint (1.3.3), I is a second-order identity tensor. 24 CHAPTER 2 GENERAL TREATMENT OF RESIDUAL STRESS IN HYPERELASTIC SPHERICAL SHELLS Section 1.3 established how residual stress would typically be expected to arise from the growth process. That growth process was characterized by decomposing the deformation gradient tensor according to (1.3.2). This made use of a reference configuration that represented a state of affairs prior to growth depicted in figure (1.4). However, for the purpose of addressing mechanical loads on the mature fully grown tissue, this is an inconvenient reference configuration. In continuum mechanics it is always possible to redefine the reference configuration. A more convenient reference configuration for this later purpose is one that describes the fully mature tissue. In making such a reference configuration change, one will now use invariants computed from a different reference configuration, and this requires a redefinition of 𝑊 so as to properly compensate. One must also account for the residual stress field. Let this residual stress field, which is due to the previous growth process, now be denoted by 𝝉. The theory of hyperelasticity with residual stress as elaborated in [85] can now be used. In the study [85], the focus was on cylindrical geometries because of their interest in tubular organs, in particular arteries. Our interest in the brain motivates a focus on spherical geometries throughout the analysis presented in Chapter 2. In this study, the effect of residual stress field 𝝉 in a hollow spherical shell on the pressure-inflation response of a Mooney-Rivlin type material is examined using methods of finite elasticity. The conception of residual stress follows that of Hoger [63, 64]. The Mooney- Rivlin form is of special interest for the sphere problem because its well-known pressure-inflation response is either monotonic or non-monotonic depending upon the shell thickness and Mooney- Rivlin parameters 𝐼1 and 𝐼2. The specific hyperelastic framework employed here follows the invariant formulation presented in [115]. That formulation is further developed and utilized in subsequent treatments for residual stress including [86], [116], and [53]. The Mooney-Rivlin form is of special interest for the sphere problem since its well-known pressure-inflation response is either monotonic or non-monotonic depending upon the shell thickness and Mooney-Rivlin parameters 𝐼1 25 and 𝐼2. The Mooney-Rivlin form of the strain-energy function is defined in terms of the invariants associated with the deformation and the residual stress fields in the reference configuration. It is assumed that the residual stress field 𝝉 resulted from the previous growth process but the residually- stressed reference configuration is not obtained by modeling the theory of finite growth explained in Section 1.3. Hence, the kinematic and mechanical quantities associated with the previous growth process such as Fe, Ce and Be can now be replaced by F, C, and B for the purpose of Chapter 2. A specific form for the residual stress component 𝜏𝑅𝑅, which is dependent only on 𝑅, is considered to define the residual stress field 𝝉. This chapter was published previously in [146], and is reproduced here with permission from Springer Nature. The present chapter is organized as follows. The general residual stress conception and its hyperelastic constitutive formulation is given in Section 2.1, followed by the specialization to spherically symmetric deformations in Section 2.2. Classes of spherically symmetric residual stress fields are introduced in Section 2.3. The integration procedure for the solution of the boundary value problem is given in Section 2.4. This leads to the consideration of specific integrals associated with the residual stress field, whose properties are examined in Section 2.5. This property determination permits a rather straight forward assessment of how the residual stress affects the pressure-inflation relation. Among the findings of interest is that particular residual stress fields can cause a monotonic inflation graph to become non-monotonic, and vice versa. Key limitations of our study are discussed in Section 2.7, some of which are associated with the elusive nature of the meaning of residual stress, along with an indication of how some of the analytical findings here could be further exploited. 2.1 Residual stress framework and constitutive modeling A body that inhabits locations X in a reference configuration B0 is considered. A residual stress is present within the body in this reference configuration, where it is given by the the tensor 𝝉. The reason for the existence of the residual stress is not important – the nature of the theory permits the consideration of residual stress fields independent of their cause. Whatever this cause, the residual stress field is subject to the classical balance of torques and forces in the usual way. Intrinsic couple 26 stresses are not considered whereupon the torque balance gives that 𝝉 is symmetric (𝝉T = 𝝉). Body forces are regarded as negligible, whereupon force balance gives that Div 𝝉 = 0, (2.1.1) where Div is the divergence operation with respect to reference locations X. The residual stress is present in the absence of surface tractions, meaning that 𝝉N = 0 on 𝜕B0, (2.1.2) where 𝜕B0 is the external boundary of B0 and N the outward pointing unit normal vector. This follows the residual stress conception of Hoger [63]. Other notions of initial stress may incorporate sustaining surface tractions, whereas none are to be present here. A consequence of conditions (2.1.1) and (2.1.2) is that the mean value of the residual stress over B0 must vanish [64], namely ∫ B0 𝝉 𝑑𝑉 = 0. (2.1.3) For this reason, any nontrivial 𝝉 will generally cause an inhomogeneous mechanical response to subsequent deformation. The application of surface tractions will deform the body to a new configuration B. Let x denote locations in B. The deformation gradient of the mapping X → x is denoted by F in the usual way. This process will generally change the internal stress field so that it is no longer given by 𝝉. Let 𝝈 denote the associated Cauchy stress tensor. If F = I then 𝝈 is equal to 𝝉 modulo a purely reactive stress contribution associated with any internal material constraints (such as incompressibility). The Cauchy stress is also symmetric and it is subject to div𝝈 = 0, (2.1.4) where div is the divergence operation with respect to deformed locations x. The surface tractions t follow from 𝝈 and the unit surface normal n on 𝜕B via t = 𝝈n on 𝜕B. (2.1.5) 27 A hyperelastic constitutive framework is employed here which follows that described in [115]. It is given in terms of a stored energy density function 𝑊 = 𝑊 (C, 𝝉) where C = FTF is the usual right Cauchy-Green deformation tensor. This form ensures the axiom of objectivity. One could more generally allow an explicit inhomogeneous response 𝑊 = 𝑊 (C, 𝝉, X) but this more general case is not considered here. It is important to note that the spatial dependence 𝝉 = 𝝉(X) confers the implicit inhomogeneous response (mentioned above) even in the context of the 𝑊 = 𝑊 (C, 𝝉) treatment. Attention is here restricted to incompressible materials, meaning that deformations X → x must be volume preserving. This gives the standard isochoric constraint The Cauchy stress now follows from F via detF = 1. 𝝈 = F 𝜕𝑊 𝜕C FT − 𝑝I, (2.1.6) (2.1.7) where 𝑝 is the non-constitutive reactive pressure associated with (2.1.6). Because F = I implies 𝝈 = 𝝉 to within a hydrostatic pressure it follows that 𝝉 = 𝜕𝑊 𝜕C (cid:12) (cid:12) (cid:12)C=I − 𝑞I, (2.1.8) for some scalar 𝑞 and, as discussed in [115], this imposes restrictions upon 𝑊. Still following [115], the dependence of 𝑊 upon C and 𝝉 can be developed in terms of a set of nine independent invariants of C, 𝝉 and their combination. These can be taken as the two usual isochoric invariants of C alone: 𝐼1 = tr C, 𝐼2 = 1 2 [(tr C)2 − tr(C2)]; (2.1.9) (note that 𝐼3 = detC = 1), the three invariants of 𝝉 alone: 𝐼41 = tr 𝝉, 𝐼42 = 1 2 [(tr 𝝉)2 − tr(𝝉2)], 𝐼43 = det 𝝉; (2.1.10) and four independent invariants of both C and 𝝉 in combination: 𝐼5 = tr(𝝉C), 𝐼6 = tr(𝝉C2), 𝐼7 = tr(𝝉2C), 𝐼8 = tr(𝝉2C2). (2.1.11) 28 Here the specific invariant sets employed in [85] are considered. The three invariants {𝐼41, 𝐼42, 𝐼43} are collectively denoted by 𝐼4. They are singled out this way because, unlike the other invariants, they do not contribute derivatives to 𝝈 via (2.1.7) by a chain rule calculation on 𝜕𝑊/𝜕C. Letting 𝑊𝑖 = 𝜕𝑊/𝜕𝐼𝑖 for 𝑖 = 1, 2, 5, 6, 7, 8 such a chain rule calculation gives 𝝈 =2W1B + 2W2(𝐼1B − B2) + 2W5𝚺 + 2W6(𝚺B + B𝚺) + 2W7𝚺B−1𝚺 + 2W8(𝚺B−1𝚺B + B𝚺B−1𝚺) − 𝑝I, (2.1.12) where, following [85], 𝚺 ≡ F𝝉FT, and B = FFT is the left Cauchy-Green tensor. Note that (2.1.12) retrieves the classical (no-residual-stress) hyperelastic result 𝝈 = 2W1B + 2W2(𝐼1B − B2) − 𝑝I when 𝝉 vanishes. A simple model form for 𝑊 to consider is one that admits the decomposition 𝑊 = 𝑊𝑜 (𝐼1, 𝐼2) + 𝑊𝜏 (𝐼5, 𝐼6, 𝐼7, 𝐼8). By taking 𝑊𝑜 to be some well vetted hyperelastic model from the conventional isotropic theory, it permits a consideration of the effect of residual stress upon standard results. Forms for 𝑊𝜏 that have received consideration (see [85]) include 𝑊𝜏 = 1 2 (𝐼5 − tr(𝝉)) and 𝑊𝜏 = 1 4 (𝐼6 − tr(𝝉)). (2.1.13) Note that 𝐼5, 𝐼6, tr(𝝉) and 𝑊𝜏 all have the same physical units, so that 1 2 and 1 numbers. The particular values 1 4 in (2.1.13) are essential as they enable consistency with 2 and 1 4 in (2.1.13) are pure (2.1.8). In what follows a Mooney-Rivlin form is considered for 𝑊𝑜, namely 𝑊𝑜 = 1 2 𝜇 (𝜅(𝐼1 − 3) + (1 − 𝜅) (𝐼2 − 3)) . (2.1.14) Here 𝜇, with units of stress, is the shear modulus. The pure number 𝜅 is taken to be in the range 0 ≤ 𝜅 ≤ 1 in order for the conventional theory to be consistent with the Baker-Ericksen inequalities. The value 𝜅 = 1 gives the neo-Hookean special case. Values for 𝜅 a bit below one are often regarded as providing a suitable first correction to the neo-Hookean theory [107]. Expressing 𝑊 = 𝑊𝑜 + 𝑊𝜏 by combining (2.1.14) with the first of (2.1.13) gives 𝑊 = 1 2 𝜇 (𝜅(𝐼1 − 3) + (1 − 𝜅) (𝐼2 − 3)) + 1 2 (𝐼5 − tr(𝝉)). (2.1.15) 29 This form for 𝑊, under the neo-Hookean specialization (𝜅 = 1), has been utilized for the analysis of residually stressed cylinders in [86] and [85]. Also, the neo-Hookean case of (2.1.15) where 𝑊𝜏 is allowed to be augmented with a quadratic term proportional to (𝐼5 − tr(𝝉))2 has been considered for the examination of acoustic waves [115] and Rayleigh waves [116] in residually stressed materials. Restricting attention to (2.1.15) it now follows from (2.1.12) that 𝝈 = 𝜇 (cid:16) (𝜅 + (1 − 𝜅)𝐼1)B − (1 − 𝜅)B2(cid:17) + F𝝉FT − 𝑝I, (2.1.16) where it is observed that F = I and 𝑝 = (2 − 𝜅)𝜇 indeed makes 𝝈 = 𝝉. 2.2 Kinematics for a radial inflation of a residually-stressed spherical shell A finite thickness spherical shell is considered in spherical coordinates (𝑅, Θ, Φ), occupies 𝐴 ≤ 𝑅 ≤ 𝐵, 0 ≤ Θ < 𝜋, 0 ≤ Φ ≤ 𝜋 (2.2.1) in the reference configuration B0. Here A and B are inner and outer radii prior to any deformation. Spherical inflation, which is one of the universal deformations in the conventional isotropic, incompressible hyperelastic theory, is then described using spherical coordinates (𝑟, 𝜃, 𝜙) in the deformed configuration B, as 𝑟 = 𝑟 (𝑅), 𝜃 = Θ, 𝜙 = Φ. (2.2.2) The deformation gradient is given by F = 𝜆𝑟e𝑟 ⊗ E𝑅 + 𝜆𝜃e𝜃 ⊗ EΘ + 𝜆𝜙e𝜙 ⊗ EΦ, (2.2.3) where {E𝑅, EΘ, EΦ} and {e𝑟, e𝜃, e𝜙} are the unit basis vectors in the reference and deformed configuration respectively, and 𝜆𝑟 = 𝑑𝑟 𝑑𝑅 , 𝜆𝜃 = 𝜆𝜙 = 𝑟 𝑅 (2.2.4) are the principal stretches in the corresponding directions. The associated right and left Cauchy- Green deformation tensors are C = 𝜆2 𝑟 E𝑅 ⊗ E𝑅 + 𝜆2 𝜃 (EΘ ⊗ EΘ + EΦ ⊗ EΦ), B = 𝜆2 𝑟 e𝑟 ⊗ e𝑟 + 𝜆2 𝜃 (e𝜃 ⊗ e𝜃 + e𝜙 ⊗ e𝜙). 30 (2.2.5) The constraint (2.1.6) determines 𝑟 (𝑅) to within a single constant parameter value, which can be taken as the radial value 𝑎 of the deformed inner radius. Then 𝑟 (𝑅) = (𝑅3 − 𝐴3 + 𝑎3)1/3, (2.2.6) whereupon the deformed sphere occupies 𝑎 ≤ 𝑟 ≤ 𝑏 = (𝐵3 − 𝐴3 + 𝑎3)1/3. Consequently the relation between 𝐴 and 𝑎 completely characterizes the amount of inflation. This also makes it convenient to emphasize the azimuthal stretch via the notation 𝜆 ≡ 𝜆𝜃 = 𝜆𝜙 = 𝑟 𝑅 ⇒ 𝜆𝑟 = 𝜆−2 = 𝑅2 𝑟 2 . (2.2.7) The reference configuration is subject to the no surface traction condition (2.1.2), and the inflation is induced by increasing the internal pressure by an amount Δ𝑃 > 0. Such loading is consistent with a radial deformation (2.2.2) in the conventional theory where there is no residual stress. In order to remain consistent with spherical inflation, consider residual stress fields of the symmetric form 𝝉 = 𝜏𝑅𝑅E𝑅 ⊗ E𝑅 + 𝜏ΘΘ(EΘ ⊗ EΘ +EΦ ⊗ EΦ) with 𝜏𝑅𝑅 = 𝜏𝑅𝑅 (𝑅) and 𝜏ΘΘ = 𝜏ΘΘ(𝑅). Then (2.1.1) and (2.1.2) respectively give and 𝑑𝜏𝑅𝑅 𝑑𝑅 + 2 𝑅 (𝜏𝑅𝑅 − 𝜏ΘΘ) = 0, 𝜏𝑅𝑅 ( 𝐴) = 0, 𝜏𝑅𝑅 (𝐵) = 0. (2.2.8) (2.2.9) A general result that follows from (2.2.8) and (2.2.9) is that ∫ 𝐵 𝐴 (cid:0)(𝑚 − 1)𝜏𝑅𝑅 + 2 𝜏ΘΘ (cid:1) 𝑅𝑚 𝑑𝑅 = 0, (2.2.10) where 𝑚 is arbitrary. The cases 𝑚 = 1 and 𝑚 = 2 then give ∫ 𝐵 𝐴 𝜏ΘΘ 𝑅 𝑑𝑅 = 0, ∫ 𝐵 𝐴 (cid:0)𝜏𝑅𝑅 + 2 𝜏ΘΘ (cid:1) 𝑅2 𝑑𝑅 = 0. (2.2.11) The first of these establishes that a nontrivial 𝜏𝜃𝜃 field must take on both tensile and compressive values. The second of these is what (2.1.3) reduces to for the problem under consideration here. 31 The Cauchy stress 𝝈 takes the form 𝝈 = 𝜎𝑟𝑟e𝑟 ⊗ e𝑟 + 𝜎𝜃𝜃 (e𝜃 ⊗ e𝜃 + e𝜙 ⊗ e𝜙) with 𝜎𝑟𝑟 = 𝜎𝑟𝑟 (𝑟) and 𝜎𝜃𝜃 = 𝜎𝜃𝜃 (𝑟). In particular, it follows from (2.1.16), (2.2.3) and (2.2.7) that 𝜎𝑟𝑟 = 𝜇 (cid:16) 2(1 − 𝜅)𝜆−2 + 𝜅𝜆−4(cid:17) + 𝜆−4𝜏𝑅𝑅 − 𝑝, and 𝜎𝜃𝜃 = 𝜇 (cid:16) (1 − 𝜅)𝜆4 + 𝜅𝜆2 + (1 − 𝜅)𝜆−2(cid:17) + 𝜆2𝜏ΘΘ − 𝑝. (2.2.12) (2.2.13) The equilibrium equation (2.1.4) gives 𝑝 = 𝑝(𝑟) along with the one nontrivial requirement 𝜕𝜎𝑟𝑟 𝜕𝑟 + 2 𝑟 (𝜎𝑟𝑟 − 𝜎𝜃𝜃) = 0. (2.2.14) In terms of the applied pressures, the surface traction condition (2.1.5) is simply 𝜎𝑟𝑟 (𝑎) = −Δ𝑃, 𝜎𝑟𝑟 (𝑏) = 0. (2.2.15) 2.3 The residual stress fields Given a 𝜏𝑅𝑅 field, the associated 𝜏ΘΘ follows from (2.2.8). Probably the simplest nontrivial 𝜏𝑅𝑅 field that is consistent with (2.2.9) is a parabolic profile proportional to (𝑅 − 𝐴) (𝑅 − 𝐵). It is expressed as 𝜏𝑅𝑅 (𝑅) = 𝛼2 𝐴2 (𝑅 − 𝐴) (𝑅 − 𝐵) = 𝛼2 𝐴2 (𝑅2 − ( 𝐴 + 𝐵)𝑅 + 𝐴𝐵) (2.3.1) where 𝛼2, with units of stress, sets the field strength. The analogous radial component field was considered in [85] for the case of a cylindrical geometry. For the spherical geometry under consideration here, using (2.2.8) with the 𝜏𝑅𝑅 field given by (2.3.1) yields 𝜏ΘΘ(𝑅) = (cid:16) 𝛼2 𝐴2 2𝑅2 − 3 2 ( 𝐴 + 𝐵)𝑅 + 𝐴𝐵 (cid:17) . (2.3.2) The subscript 2 in 𝛼2 is indicative of the highest power of 𝑅 in these expressions. We wish to consider a wider variety of residual stress fields (𝜏𝑅𝑅 (𝑅), 𝜏ΘΘ(𝑅)) beyond that given by (2.3.1) and (2.3.2). For reasons that will become clear later, it will also be advantageous to consider 𝜏𝑅𝑅 fields that consist of a linear combination of terms in 𝑅0, 𝑅2, 𝑅5, and more generally 𝑅2+𝑛 for 𝑛 = 0, 1, 2, . . . . Then the boundary conditions (2.2.9) motivates the consideration of 𝜏𝑅𝑅 (𝑅) = 𝛼5 𝐴5 (cid:16) 𝑅5 − (cid:17) (cid:16) 𝐵5−𝐴5 𝐵2−𝐴2 𝑅2 + 𝐴2𝐵2 (𝐵3−𝐴3) 𝐵2−𝐴2 (cid:17) , (2.3.3) 32 as well as higher order forms such as and 𝜏𝑅𝑅 (𝑅) = 𝛼8 𝐴8 (cid:16) 𝑅8 − (cid:17) (cid:16) 𝐵8−𝐴8 𝐵5−𝐴5 𝑅5 + 𝐴5𝐵5 (𝐵3−𝐴3) 𝐵5−𝐴5 (cid:17) , 𝜏𝑅𝑅 (𝑅) = (cid:16) 𝛼11 𝐴11 𝑅11 − (cid:16) 𝐵11−𝐴11 𝐵8−𝐴8 (cid:17) 𝑅8 + 𝐴8𝐵8 (𝐵3−𝐴3) 𝐵8−𝐴8 (cid:17) . (2.3.4) (2.3.5) The multiplier coefficients 𝛼2, 𝛼5, 𝛼8, 𝛼11 in (2.3.1)-(2.3.5), all with units of stress, will serve as a useful marker for distinguishing between these different fields in the various examples that follow. The 𝜏ΘΘ corresponding to (2.3.3)-(2.3.5) will again follow from (2.2.8). In making such a computation, and especially for later developments, it is useful to proceed term by term in the residual stress expressions such as (2.3.3). The overbar notation ¯𝜏 is utilized for this purpose. Thus employing (2.2.8) gives the single term result ¯𝜏𝑅𝑅 = 𝑐 𝑅𝑞, ⇒ ¯𝜏ΘΘ = 𝑐 (cid:16) 1 + 1 2 𝑞 (cid:17) 𝑅𝑞. Consequently, if 𝜏𝑅𝑅 is given by (2.3.3) then 𝜏ΘΘ(𝑅) = 𝛼5 𝐴5 (cid:16) 7 2 𝑅5 − 2 (cid:17) (cid:16) 𝐵5−𝐴5 𝐵2−𝐴2 𝑅2 + 𝐴2𝐵2 (𝐵3−𝐴3) 𝐵2−𝐴2 (cid:17) . (2.3.6) (2.3.7) The purpose of the ¯𝜏 notation is to provide a reminder that a single term within the overall multi-term expression is unlikely to be consistent with (2.2.9). Thus (2.3.7) is consistent with (2.2.9) whereas (2.3.6) is not. Rather, (2.3.6) is a preliminary single term result that will be used in conjunction with other single term expressions to achieve consistency with the boundary conditions. For completeness in what follows, the 𝜏ΘΘ field associated with (2.3.4) now also follows as 𝜏ΘΘ(𝑅) = (cid:16) 𝛼8 𝐴8 5𝑅8 − 7 2 (cid:17) (cid:16) 𝐵8−𝐴8 𝐵5−𝐴5 𝑅5 + 𝐴5𝐵5 (𝐵3−𝐴3) 𝐵5−𝐴5 (cid:17) , and the 𝜏ΘΘ field associated with (2.3.5) follows as 𝜏ΘΘ(𝑅) = 𝛼11 𝐴11 (cid:16) 13 2 𝑅11 − 5 (cid:16) 𝐵11−𝐴11 𝐵8−𝐴8 (cid:17) 𝑅8 + 𝐴8𝐵8 (𝐵3−𝐴3) 𝐵8−𝐴8 (cid:17) . (2.3.8) (2.3.9) Figure 2.1 shows the effect of thickness on the residual stress field pair given by (2.3.1) and (2.3.2). The 2-D counterpart to this figure is Figure 1 of [85]. As in the 2-D case, the 𝜏ΘΘ fields 33 Figure 2.1 Residual stress fields (2.3.1) for 𝜏𝑅𝑅 (red) and (2.3.2) for 𝜏ΘΘ (blue) showing the effect of shell thickness. Here 𝛼2 = 1, 𝐴 = 1 and: (a) B=1.5, (b) B=2, (c) B=3, and (d) B=4. transition from monotone to non-monotone as the thickness increases. The minimum develops sooner in the present 3-D case since, unlike the 2-D case, it is readily apparent here for the case of 𝐴 = 1, 𝐵 = 2. Note also the qualitative similarity of these fields, which are independent of constitutive law, to the residual stress field displayed in Figure 3 of [73] that is associated with the eversion of a spherical shell composed of a particular Mooney-Rivlin material. For a common shell thickness, Figure 2.2 displays the four different stress field pairs (𝜏𝑅𝑅, 𝜏ΘΘ) as 𝜏𝑅𝑅 is given by the four alternatives: (2.3.1), (2.3.3), (2.3.4) and (2.3.5). In particular, 𝜏𝑅𝑅 is of a single sign in all four cases. In contrast, 𝜏ΘΘ takes on both positive and negative values, a result which follows from the first of (2.2.11). Note also that the fields increasingly localize to the outer radius 𝑅 = 𝐵 as the order (highest power of 𝑅) increases in the residual stress fields. This 34 Figure 2.2 Residual stress fields 𝜏𝑅𝑅 in red and 𝜏ΘΘ in blue showing how the profiles change as higher order radial terms are employed. Here 𝐴 = 1, 𝐵 = 2 and 𝛼𝑖 = 1 in all four cases. Panel (a) is for (2.3.1) and (2.3.2); panel (b) is for (2.3.3) and (2.3.7); panel (c) is for (2.3.4) and (2.3.8); panel (d) is for (2.3.5) and (2.3.9). tendency is highlighted in Figure 2.3 which redisplays the stress fields from Figure 2.2, but now with 𝛼𝑖 chosen so that 𝜏ΘΘ(𝐵) = 1 for all four alternatives. 2.4 Pressure-inflation relations For a given residually stressed spherical shell, one seeks to determine the relation between applied load (embodied in 𝑃𝑎 and 𝑃𝑏) and the amount of inflation (the value of 𝑎 in relation to 𝐴). For this purpose (2.2.14) admits to the integration 2(𝜎𝜃𝜃 − 𝜎𝑟𝑟) 𝑑𝑟 𝑟 . (2.4.1) 𝜎𝑟𝑟 (𝑟) − 𝜎𝑟𝑟 (𝑎) = ∫ 𝑟 𝑎 35 Figure 2.3 Stress fields from Figure 2.2 with 𝜏𝑅𝑅 in the left panel and 𝜏ΘΘ in the right panel, now normalized by taking 𝛼2 = 1, 𝛼5 = 3/116, 𝛼8 = 31/11344, 𝛼11 = 225/776192. This makes 𝜏ΘΘ(𝐵) = 1. Evaluation of this at 𝑟 = 𝑏 using the boundary conditions (2.2.15) gives Δ𝑃 = ∫ 𝑏 𝑎 2(𝜎𝜃𝜃 − 𝜎𝑟𝑟) 𝑑𝑟 𝑟 . (2.4.2) Substituting from (2.2.12) and (2.2.13) into the above expression eliminates the hydrostatic pressure 𝑝 from the treatment. It is further convenient to then decompose the right side of (2.4.2) as ∫ 𝑏 𝑎 2(𝜎𝜃𝜃 − 𝜎𝑟𝑟) 𝑑𝑟 𝑟 = 𝑃𝑜 + 𝑃𝜏, (2.4.3) where 𝑃𝑜 and 𝑃𝜏 gather together the contributions arising from 𝑊𝑜 and 𝑊𝜏, respectively. Conse- quently with and Δ𝑃 = 𝑃𝑜 + 𝑃𝜏 𝑃𝑜 = 2𝜇 ∫ 𝑏 (cid:18) 𝑎 𝜅(𝜆2 − 𝜆−4) + (1 − 𝜅) (𝜆4 − 𝜆−2) (cid:19) 𝑑𝑟 𝑟 , 𝑃𝜏 = 2 ∫ 𝑏 (cid:18) 𝑎 𝜆2𝜏ΘΘ − 𝜆−4𝜏𝑅𝑅 (cid:19) 𝑑𝑟 𝑟 . (2.4.4) (2.4.5) (2.4.6) 36 The term 𝑃𝑜 is familiar from the conventional hyperelastic theory and readily admits integration using 𝑑𝑟/𝑟 = 𝑑𝜆/(𝜆 − 𝜆4). Introduce non-dimensionalized quantities 𝑃∗ 𝑜 = 𝑃𝑜 𝜇 , 𝜂 = 𝐵 𝐴 , 𝑒 = 𝑎3 𝐴3 − 1, (2.4.7) and note that 𝜂 serves as an aspect ratio or relative thickness parameter that characterizes the initial shell geometry. The parameter 𝑒 represents the relative inflation, with 𝑒 = 0 corresponding to the reference state and 𝑒 > 0 corresponding to inflation. The integration of (2.4.5) now yields 𝑃∗ 𝑜 = 𝜅 (cid:19) − (cid:18) 5𝜂4 + 4𝑒𝜂 2(𝜂3 + 𝑒) 4 4𝑒 + 5 2(𝑒 + 1) 4 (cid:18) 1 + 2𝑒 (𝑒 + 1) 2 𝑜 vanishes when 𝑒 = 0, which is consistent with zero pressure in the 𝜂3 + 2𝑒 𝜂(𝜂3 + 𝑒) 2 + (1 − 𝜅) (2.4.8) − (cid:19) . 3 3 3 3 It should be noted that 𝑃∗ reference (undeformed) state. The inflation relation for non-residually stressed hyperelastic spherical shells has been the subject of longstanding study. Hence, a quick recap of the pressure-inflation response for Mooney- Rivlin spherical shells is given in the rest of the chapter, as given by (2.4.8). In other words, this recap recalls well known results because it is for the case in which there is no residual stress. The associated pressure-inflation response graphs, with Δ𝑃 on the ordinate and either 𝑒 or 𝑎/𝐴 on the abscissa, all exhibit an initial inflation, meaning that they all start from from Δ𝑃 = 0 when 𝑒 = 0 or 𝑎/𝐴 = 1, and all such graphs initially increase with 𝑒. However as 𝑒 continues to increase, one of three behavior types is possible, depending on the values of 𝜅 and 𝜂. Type (a) behavior is one of continued monotone Δ𝑃 increase as 𝑒 → ∞. Type (b) behavior involves increase to a local maximum followed by monotone decrease with Δ𝑃 → 0 in an asymptotic fashion as 𝑒 → ∞. Type (c) behavior involves three separate intervals of Δ𝑃 response; the first interval involves increase to a local maximum, followed by decrease to a local minimum (with Δ𝑃 > 0), and then concluding with monotone graphical increase as 𝑒 → ∞. Type (b) behavior is restricted to the neo-Hookean case (𝜅 = 1) and this behavior type occurs for all aspect ratios 𝜂 when the material is neo-Hookean. Type (a) behavior will occur if the Mooney- Rivlin parameter 𝜅 < 0.823 (an exact root expression involving cubics will provide additional precision [26]), and this is again true for all aspect ratios 𝜂. For 0.823 < 𝜅 < 1, either type 37 (a) or type (c) behavior is possible depending on the aspect ratio 𝜂. If the shell is sufficiently thick then type (a) behavior occurs. If it is sufficiently thin then type (c) behavior occurs. The relation between 𝜂 and 𝜅 that delineates between these two behaviors can be characterized exactly and elegantly, not only for the Mooney-Rivlin material, but also for more general incompressible isotropic hyperelastic materials (in which case additional behavior types are not ruled out [26]). For our purposes in what follows, we shall continue to provide certain demonstrations using the Mooney-Rivlin base response for a spherical shell with 𝜂 = 𝐵/𝐴 = 2. Then, using the procedure as described in [149], one finds that the transition between type (a) and type (c) behavior occurs for 𝜂 = 2 when 𝜅 (cid:27) 0.847526. This is illustrated in Figure 2.4 where, for 𝜂 = 2, one finds that: 𝜅 = 0.8 gives type (a) behavior, 𝜅 = 0.9 gives type (c) behavior, and 𝜅 = 1.0 gives type (b) behavior. Figure 2.4 Mooney-Rivlin material response for a spherical shell with 𝜂 = 𝐵/𝐴 = 2 showing the effect of 𝜅. The neo-Hookean case 𝜅 = 1 gives type (b) behavior. Type (a) response occurs for 0 ≤ 𝜅 < 0.848 as shown by the 𝜅 = 0.8 response curve. Type (c) response occurs for 0.848 < 𝜅 < 1 as shown by the 𝜅 = 0.9 response curve. Left panel shows pressure as a function of 𝑒 and right panel shows pressure as a function of 𝑎/𝐴. 2.5 The residual stress integrals The residual stress effect is governed by the integral for 𝑃𝜏 in (2.4.6). Because the residual stress fields under consideration occur as linear combinations of powers of 𝑅, it is useful to first 38 observe that ¯𝜏𝑅𝑅 = 𝑐 𝑅𝑞 ⇒ 𝜆2 ¯𝜏ΘΘ − 𝜆−4 ¯𝜏𝑅𝑅 = 𝑐 (cid:16) 𝜆2(1 + 1 2 𝑞) − 𝜆−4(cid:17) 𝑅𝑞. This motivates the definition 𝑃 ¯𝜏|𝑞 ≡ ∫ 𝑏 𝑎 (cid:16) (2 + 𝑞)𝜆2 − 2𝜆−4(cid:17) 𝑅𝑞 𝑑𝑟 𝑟 . (2.5.1) The integral for 𝑃𝜏 in (2.4.6) will now be expressed in terms of these 𝑃 ¯𝜏|𝑞. For example, the residual stress field (2.3.1), (2.3.2) makes 𝑃𝜏 = (cid:16) 𝛼2 𝐴2 𝑃 ¯𝜏|2 − ( 𝐴 + 𝐵)𝑃 ¯𝜏|1 + 𝐴𝐵 𝑃 ¯𝜏|0 (cid:17) . (2.5.2) 2.5.1 Explicit integration of 𝑃 ¯𝜏|𝑞 The explicit integration of 𝑃 ¯𝜏|𝑞 is complicated by the fact that now 𝑅 appears in the integrand, along with 𝜆 and 𝑟. This prompts the introduction of dummy integration variable 𝛽 = (𝜆3 − 1)1/3 whereupon 𝑅 = 𝑒1/3 𝐴 𝛽 , 𝜆 = (𝛽3 + 1)1/3, 𝑑𝑟 𝑟 = − 𝑑𝛽 𝛽4 + 𝛽 . This motivates the additional non-dimensionalized quantity 𝑃∗ ¯𝜏|𝑞 = 𝑃 ¯𝜏|𝑞 𝐴𝑞 . One now finds that 𝑃∗ ¯𝜏|𝑞 = −2𝑒 𝑞 3 ∫ 𝑒 1 3 /𝜂 1 3 𝑒 (cid:18) 𝑞 + 2 2 − 1 (𝛽3 + 1)2 (cid:19) 𝑑𝛽 (𝛽3 + 1) 1 3 𝛽𝑞+1 (2.5.3) (2.5.4) (2.5.5) Whether or not (2.5.5) integrates simply depends upon the value of 𝑞. Simple integration, in the sense of not requiring representation in terms of special functions, occurs for 𝑞 = 0, 2, 5, 8, . . . 2 + 3𝑛, . . . In particular, we find that 𝑃∗ ¯𝜏|0 = 5 + 4𝛽3 2(1 + 𝛽3) 4 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 39 (2.5.6) ¯𝜏|2 = 𝑒 2 𝑃∗ 3 2 − 4𝛽3 − 5𝛽6 2𝛽2(1 + 𝛽3) 4 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 ¯𝜏|5 = 𝑒 5 𝑃∗ 3 2 + 5𝛽3 + 16𝛽6 + 12𝛽9 2𝛽5(1 + 𝛽3) 4 3 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) ¯𝜏|8 = 𝑒 8 𝑃∗ 3 10 + 16𝛽3 − 22𝛽6 − 132𝛽9 − 99𝛽12 10𝛽8(1 + 𝛽3) 4 3 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) and (2.5.7) (2.5.8) (2.5.9) ¯𝜏|11 = 𝑒 11 𝑃∗ 3 40 + 55𝛽3 − 28𝛽6 + 126𝛽9 + 756𝛽12 + 567𝛽15 40𝛽11(1 + 𝛽3) 4 3 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (2.5.10) However other values of 𝑞 do not lead to such simple expressions. For example 𝑞 = 1 gives ¯𝜏|1 =𝑒 1 𝑃∗ 3 (cid:16) 𝛽−1 2𝐹1 (cid:2) − 1 3 , 7 3 , 2 3 , −𝛽3(cid:3) − 3𝛽2 2𝐹1 (cid:2) 2 3 , 7 3 , 5 3 , −𝛽3(cid:3) 𝛽5 − 3 5 2𝐹1 (cid:2) 5 3 , 7 3 , 8 3 , −𝛽3(cid:3) (cid:17) 𝛽=𝑒 1 3 /𝜂 𝛽=𝑒 1 3 (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (2.5.11) where 2𝐹1 [·, ·, ·, ·] is the hypergeometric function (see, e.g., [2]) that is defined by the series 2𝐹1(𝑎, 𝑏, 𝑐, 𝑧) = Γ (𝑐) Γ (𝑎) Γ (𝑏) ∞ ∑︁ 𝑛=0 Γ (𝑎 + 𝑛) Γ (𝑏 + 𝑛) Γ (𝑐 + 𝑛) 𝑧𝑛 𝑛! . (2.5.12) Here Γ(·) is the classical gamma function, which has many formal expressions [95], including Γ(𝜉) = ∫ ∞ 0 𝑡𝜉−1 𝑒−𝑡 𝑑𝑡. (2.5.13) All of the 𝑃∗ positive values. Beginning from the initial value 𝑃∗ ¯𝜏|𝑞 as functions of 𝑒 for 𝑒 > 0 as given by (2.5.6) - (2.5.11) are found to take on = 𝜂𝑞 − 1, the graphs are found to increase (cid:12) (cid:12) (cid:12)𝑒=0 ¯𝜏|𝑞 40 to a maximum and to then subsequently decrease toward zero as 𝑒 tends to infinity. This behavior is established in the next subsection. 2.5.2 Behavior of 𝑃∗ ¯𝜏|𝑞 for small and large 𝑒 It is useful to obtain simplified expressions for both the small 𝑒 and large 𝑒 behavior of 𝑃∗ ¯𝜏|𝑞 as functions of the parameters 𝜂 and 𝑞. This is especially true for those 𝑞, such as 𝑞 = 1, that do not give straight forward expressions for 𝑃∗ ¯𝜏|𝑞. To obtain the small 𝑒 response we rewrite (2.5.5) as with 𝑃∗ ¯𝜏|𝑞 = 𝑒 𝑞 3 ∫ 𝑒 1 3 /𝜂 1 3 𝑒 𝑔(𝑞, 𝛽) 𝑑𝛽 𝛽𝑞+1 𝑔(𝑞, 𝛽) = 1 (𝛽3 + 1) 1 3 (cid:18) −𝑞 − 2 + (cid:19) . 2 (𝛽3 + 1)2 (2.5.14) (2.5.15) Because 𝑔(𝑞, 𝛽) has the small 𝛽 expansion 𝑔(𝑞, 𝛽) = −𝑞 + 𝑎1𝛽3 + 𝑎2𝛽6 + 𝑎3𝛽9 + 𝑂 (𝛽12) (2.5.16) with 𝑎1 = −4 + 1 3 𝑞, it follows from (2.5.14) that 𝑎2 = 22 3 − 2 9 𝑞, 𝑎3 = 46 9 + 14 81 𝑞, ¯𝜏|𝑞 = 𝑒𝑞/3 𝑃∗ (cid:18) 𝛽−𝑞 + 𝑎1 3 − 𝑞 𝛽3−𝑞 + 𝑎2 6 − 𝑞 𝛽6−𝑞 + 𝑎3 9 − 𝑞 𝛽9−𝑞 + 𝑂 (𝛽12−𝑞) (cid:19) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) 𝛽=𝑒 1 3 /𝜂 , 𝛽=𝑒 1 3 (2.5.17) provided 𝑞 ≠ 0, 3, 6, 9. The cases 𝑞 = 0, 3, 6, 9 . . . can be considered separately, with 𝑞 = 3, 6, 9 . . . each yielding up an alternative log term in the above expression. Finally, evaluating (2.5.17) between the upper and lower limits now shows that the small 𝑒 expansion for 𝑃∗ ¯𝜏|𝑞 begins as ¯𝜏|𝑞 = 𝜂𝑞 − 1 − 𝑃∗ (cid:19) (cid:18) 12 − 𝑞 3 − 𝑞 1 3 (𝜂𝑞−3 − 1)𝑒 + 𝑂 (𝑒2), (2.5.18) provided 𝑞 ≠ 3 (with 𝑞 = 3 easily treated as its own special case). 41 For large 𝑒 consider the large 𝛽 behavior of 𝑔(𝑞, 𝛽) using (2.5.15), namely 𝑔(𝑞, 𝛽) ∼ −(𝑞 + 2)/𝛽 as 𝛽 → ∞. This gives 𝑃∗ ¯𝜏|𝑞 ∼ −(𝑞 + 2)𝑒 𝑞 3 ∫ 𝑒 1 3 /𝜂 1 3 𝑒 𝑑𝛽 𝛽𝑞+2 = (cid:19) (cid:18) 𝑞 + 2 𝑞 + 1 (𝜂𝑞+1 − 1)𝑒−1/3, as 𝑒 → ∞. (2.5.19) The first three panels of Figure 2.5 – which consider separately 𝑞 = 0, 1, 2 – provides graphical verification that the small 𝑒 and large 𝑒 results give good approximations to the exact integrals on an appropriate range. The final panel of this figure puts all of the graphs for (2.5.6) - (2.5.11) on a common plot, so as to also show results for the remaining cases 𝑞 = 5, 8, 11. 2.5.3 Consequences for 𝑃𝜏 The 𝑃𝜏 integral for each of the previously exhibited 𝜏𝑅𝑅 fields (2.3.1), (2.3.3), (2.3.4), (2.3.5) are expressed in terms of the nondimensionalized 𝑃∗ ¯𝜏|𝑞 as follows: (2.3.1) ⇒ 𝑃𝜏 = 𝛼2 (cid:16) ¯𝜏|2 − (1 + 𝜂) 𝑃∗ 𝑃∗ ¯𝜏|1 + 𝜂 𝑃∗ ¯𝜏|0 (cid:17) , (2.3.3) ⇒ 𝑃𝜏 = 𝛼5 (2.3.4) ⇒ 𝑃𝜏 = 𝛼8 (cid:16) (cid:16) 𝑃∗ ¯𝜏|5 − (cid:17) (cid:16) 𝜂5−1 𝜂2−1 𝑃∗ ¯𝜏|2 + (cid:16) 𝜂2 (𝜂3−1) 𝜂2−1 (cid:17) (cid:17) , 𝑃∗ ¯𝜏|0 𝑃∗ ¯𝜏|8 − (cid:17) (cid:16) 𝜂8−1 𝜂5−1 𝑃∗ ¯𝜏|5 + (cid:16) 𝜂5 (𝜂3−1) 𝜂5−1 (cid:17) (cid:17) , 𝑃∗ ¯𝜏|0 (2.3.5) ⇒ 𝑃𝜏 = 𝛼11 (cid:16) 𝑃∗ ¯𝜏|11 − (cid:17) (cid:16) 𝜂11−1 𝜂8−1 𝑃∗ ¯𝜏|8 + (cid:16) 𝜂8 (𝜂3−1) 𝜂8−1 (cid:17) (cid:17) . 𝑃∗ ¯𝜏|0 (2.5.20) (2.5.21) (2.5.22) (2.5.23) The coefficients 𝛼𝑖 in (2.5.20) - (2.5.23) have units of stress. Here it is important to note (see, e.g., Figure 2.2) that setting 𝛼𝑖 = 1 corresponds to very different residual stress field magnitudes for the four cases 𝑖 = 1, . . . 4. In order to fairly contrast the different fields, we observe that 𝜏ΘΘ(𝐵) is a dominant residual stress value in Figure 2.2. This motivates introducing the non-dimensional 𝛼∗ 𝑖 = 𝜏ΘΘ(𝐵)/𝜇. Consequently, using (2.3.2), (2.3.7), (2.3.8), (2.3.9), one finds that 𝛼2 = 2 𝜂 (𝜂 − 1) 𝜇, 𝛼∗ 2 𝛼5 = 2(𝜂2 − 1) 𝜂2 (3𝜂5 − 5𝜂3 + 2) 𝜇, 𝛼∗ 5 (2.5.24) (2.5.25) 42 Figure 2.5 The solid curves in panels (a) - (c) show the 𝑃∗ ¯𝜏|𝑞 response graphs as a function of 𝑒 for 𝜂 = 2. Panel (a) is for 𝑞 = 0 as given by (2.5.6). Panel (b) is for 𝑞 = 1 as given by (2.5.11). Panel (c) is for 𝑞 = 2 as given by (2.5.7). The dashed curves in each of these panels depicts the associated small 𝑒 (red) and large 𝑒 (blue) response as given by (2.5.18) and (2.5.19), respectively. The final panel (d) shows all six graphs 𝑃∗ ¯𝜏|𝑞 for 𝑞 = 0, 1, 2, 5, 8, 11 on a common plot. Because the initial = 𝜂𝑞 − 1 = 2𝑞 − 1, only the cases 𝑞 = 8 and 𝑞 = 11 are well resolved in the last panel, value 𝑃∗ with the other cases crowded near the 𝑒-axis. (cid:12) (cid:12) (cid:12)𝑒=0 ¯𝜏|𝑞 43 0500.51 𝛼8 = 2(𝜂5 − 1) 𝜂5 (3𝜂8 − 8𝜂3 + 5) 𝜇, 𝛼∗ 8 𝛼11 = 2(𝜂8 − 1) 𝜂8 (3𝜂11 − 11𝜂3 + 8) 𝛼∗ 11 𝜇. (2.5.26) The final quantities to nondimensionalize are the 𝑃𝜏 in (2.5.20) - (2.5.23). The previous nondimensionalization 𝑃∗ 𝑜 = 𝑃𝑜/𝜇 from the first of (2.4.7) motivates 𝑃∗ 𝜏 = 𝑃𝜏/𝜇 whereupon (2.4.4) is expressed in the nondimensionalized form Δ𝑃/𝜇 = 𝑃∗ 𝑜 + 𝑃∗ 𝜏, (2.5.27) (2.5.28) where 𝑃∗ 𝑜 continues to be given by (2.4.8). More importantly, 𝑃∗ 𝜏 in (2.5.28) now follows from (2.5.20) - (2.5.23) in terms of 𝛼∗ 𝑖 . For example, (2.3.1) ⇒ 𝑃∗ 𝜏 = 2 𝛼∗ 2 𝜂(𝜂 − 1) (cid:16) ¯𝜏|2 − (1 + 𝜂) 𝑃∗ 𝑃∗ ¯𝜏|1 + 𝜂 𝑃∗ ¯𝜏|0 (cid:17) , (2.5.29) (2.3.3) ⇒ 𝑃∗ 𝜏 = 2(𝜂2 − 1) 𝛼∗ 5 𝜂2(3𝜂5 − 5𝜂3 + 2) (cid:16) 𝑃∗ ¯𝜏|5 − (cid:17) (cid:16) 𝜂5−1 𝜂2−1 𝑃∗ ¯𝜏|2 + (cid:16) 𝜂2 (𝜂3−1) 𝜂2−1 (cid:17) (cid:17) . 𝑃∗ ¯𝜏|0 (2.5.30) Similar style expressions, albeit more complicated ones, follow for 𝑃∗ 𝜏 upon using (2.5.26) in the remaining two cases of (2.5.22) and (2.5.23). Together this gives four expressions for 𝑃∗ 𝜏: one each for (2.3.1), (2.3.3), (2.3.4) and (2.3.5). Each of these is a function of 𝑒 on 𝑒 > 0 as determined by the parameters 𝜂 > 1 and 𝛼∗ 𝑖 . In each case, and for all parameter values 𝜂 and 𝛼∗ 𝑖 , one may verify that (2.5.18) gives that 𝑃∗ 𝜏 = 0 when 𝑒 = 0. The advantage of introducing the unit-free coefficient 𝛼∗ 𝑖 is that the various alternative residual stress fields are now expressed in terms of the readily interpretable scaling ratio 𝜏ΘΘ(𝐵)/𝜇. For 𝑖 > 0 we find that each of the 𝑃∗ 𝛼∗ 𝜏 decreases with 𝑒 from its initial zero value to a local minimum before again increasing so as to asymptotically give 𝑃∗ 𝜏 → 0 as 𝑒 → ∞. The first panel of Figure 2.6 shows this behavior for (2.3.1) and (2.3.4) on the interval 0 ≤ 𝑒 ≤ 25. The panel shows how this behavior is also reliably predicted upon using the small and large 𝑒 results for the individual parts 𝑃∗ ¯𝜏|𝑞 that make up each 𝑃∗ 𝜏. The second panel of Figure 2.6 confirms that these trends hold 44 𝜏 as a function of 𝑒 for the four cases (2.3.1), (2.3.3), (2.3.4) and (2.3.5). Figure 2.6 Behavior of 𝑃∗ All graphs are for 𝜂 = 2 and 𝛼∗ 𝑖 = 1. Panel (a) shows the case (2.3.1) in purple and the case (2.3.4) in green. Small and large 𝑒 approximations are shown in the same color using dashed lines. Panel (b) supplies the remaining graphs for (2.3.3) and (2.3.5). The general behavior for 𝛼∗ 𝑖 > 0 is monotone decrease to a local minimum followed by asymptotic increase back toward the 𝑒-axis. Changing the sign of 𝛼∗ 𝑖 < 0 reflects all curves about the horizontal 𝑒-axis. 𝑖 so as to make 𝛼∗ also for the remaining two cases by showing 𝑃∗ 𝜏 on the interval 0 ≤ 𝑒 ≤ 40 for all four cases (2.3.1), (2.3.3), (2.3.4) and (2.3.5). The new scaling (2.5.24) now results in the case (2.3.1) providing the dominant affect for a common value of the dimensionless ratio 𝜏ΘΘ(𝐵)/𝜇. 2.6 Effect of residual stress on the pressure-inflation relation The vanishing of 𝑃∗ 𝜏 both at 𝑒 = 0 and as 𝑒 → ∞ indicates that the residual stress effect on the inflation behavior is most consequential for intermediate values of 𝑒 for the residual stress fields under consideration here. For example, the basic 𝜏𝑅𝑅 and 𝜏ΘΘ field pair given by (2.3.1) and (2.3.2) generates (2.5.29) for 𝑃∗ 𝜏. If 𝛼∗ 2 > 0 then this 𝑃∗ 𝜏 is negative, with an 𝜂 dependent value of 𝑒 that locates its internal minimum. For 𝜂 = 2 this internal minimum is at 𝑒 (cid:27) 3.9543. Consequently, it is at this 𝑒-value that the residual stress will most affect the inflation graph. Figure 2.7 exhibits this tendency, where for demonstration purposes the Mooney-Rivlin parameter 𝜅 = 0.75 is employed. The baseline curve for no residual stress is type (a) for 𝜅 = 0.75. As confirmed by the figure, the effect of the residual stress field (2.3.1), (2.3.2) on the inflation graph is largest near 𝑒 = 4. Positive values of 𝛼∗ 2 correspond to 𝜏𝑅𝑅 < 0 in the shell interior and these lower the inflation graphs. Conversely, interior residual stresses 𝜏𝑅𝑅 > 0 raise the graph. This reflects the same tendency for 45 the corresponding 2-D cylindrical geometry problem described in [85]. Similar results follow for the other 𝜏𝑅𝑅 fields (2.3.3), (2.3.4) and (2.3.5), although the overall effect is smaller for the same value of 𝜏ΘΘ(𝐵) because of the overall ordering of the various 𝑃∗ 𝜏 graphs as previously exhibited in Figure 2.6. Figure 2.7 Inflation graphs as determined by (2.5.28) for a 𝜂 = 𝐵/𝐴 = 2 spherical shell consisting of a hyperelastic base material described by a 𝜅 = 0.75 Mooney-Rivlin stored energy density 𝑊. The residual stress field is given by (2.3.1) and (2.3.2). The affect of the residual stress is accounted for in the theory via (2.1.15) which augments the Mooney-Rivlin energy density with the term 2 ((tr(𝝉𝐶) − tr(𝝉)). The central curve in red is the base inflation response (𝛼∗ 1 2 = 0) corresponding to no residual stress. Solid blue curves correspond progressively to 𝛼∗ 2 = 0.5, 1, 2. Dashed blue curves correspond progressively to 𝛼∗ 2 = −0.5, −1, −2. The specific nature of (2.1.15), with its separate additive contributions, accounts for the straight forward additive residual stress effect in (2.5.28) that serves to modify the base response graph (the red curve in Figure 2.7). In particular, for 𝜂 = 2 and stress field (2.3.1), (2.3.2), the residual 46 stress will continue to have its largest effect on the inflation graph near 𝑒 = 4 for other values of the Mooney-Rivlin parameter 𝜅. Indeed, the same holds true if 𝑊 in (2.1.15) is modified so as to replace the Mooney-Rivlin form with an alternative base material response. Staying with the Mooney-Rivlin base response, but considering alternative values of 𝜅, recall from the discussion at the end of section 2.4 that the value 𝜅 (cid:27) 0.847526 provides the transition between type (a) and type (c) inflation response for 𝜂 = 2 when no residual stress is present. The question thus arises as to whether residual stress can alter the type of the inflation response away from its base behavior? The answer to this question is “Yes". As shown in our final two figures, residual stress can cause a base type (a) response to become type (c), and it can also cause a base type (c) response to become type (a). Technically, such a result can be made to follow from the decreasing-increasing behavior of 𝑃∗ 𝜏 for 𝛼∗ 𝑖 > 0 and its converse increasing-decreasing behavior for 𝛼∗ 𝑖 < 0. Namely, prior to the consideration of residual stress, the transition between type (a) and type (c) behavior is characterized by an inflation graph that is locally flat at some value of 𝑒. Specifically, this 𝑒 locates an inflection point of the base response graph. By contriving to make the inflection point of the base response graph have an 𝑒-value that is suitably located with respect to the single extrema of a 𝑃∗ 𝜏 graph, one can chose 𝛼∗ 𝑖 to push the base response graph into either type (a) or type (c) behavior. More generally, starting with a base response graph that is just to one side of the (a)-to-(c) transition, namely either a definite type (a) response or a definite type (c) response, one can by continuity execute a similar push so as to transition the inflation graph to the other type of response behavior. To provide a demonstration of a type (a) base response that becomes type (c) for a specific residual stress field, consider 𝜂 = 2 and 𝜅 = 0.84. This corresponds to type (a) response in the absence of residual stress (the red curve in Figure 2.8). The residual stress field pair given by (2.3.1) and (2.3.2) with 𝛼∗ 2 = −1 is then found to give a type (c) inflation response graph (the dashed blue curve in Figure 2.8). Conversely we exhibit a type (c) base response that becomes type (a) for a specific residual stress field by taking 𝜂 = 2 and 𝜅 = 0.86. This corresponds to type (c) response in the absence of residual stress (the red curve in Figure 2.9). The residual stress field pair given by 47 Figure 2.8 Type (a) monotonic inflation response occurring for 𝜂 = 2 and 𝜅 = 0.84 in the absence of residual stress (red) becomes type (c) nonmonotonic response for the residual stress field pair given by (2.3.1) and (2.3.2) with 𝛼∗ 2 = −1 (dashed blue). (2.3.1) and (2.3.2) with 𝛼∗ 2 = 1 is then found to give a type (a) inflation response graph (the solid blue curve in Figure 2.9). 2.7 Additional remarks The specific additive form 𝑊 = 𝑊𝑜 +𝑊𝜏 that is used in this study is motivated by, and in keeping with, prototypical forms that are utilized in previous works. These additive forms are now beginning to see usage in computational treatments (e.g., [110, 41]) that address fundamental questions of bifurcation, imperfection sensitivity and symmetry breaking deformations, all of which are issues that are beyond the scope of this article. In the context of 𝑊 = 𝑊𝑜 + 𝑊𝜏 we have made use of 𝑊𝜏 given by the first of (2.1.13). For spherical symmetry this gave rise to integrals (2.4.6) which had readily determined properties for certain residual stress field expressions. Specifically, a 𝜏𝑅𝑅 field 48 Figure 2.9 Type (c) nonmonotonic inflation response occurring for 𝜂 = 2 and 𝜅 = 0.86 in the absence of residual stress (red) becomes type (a) monotonic response for the residual stress field pair given by (2.3.1) and (2.3.2) with 𝛼∗ 2 = 1 (blue). that is a linear combination of terms in 𝑅𝑞 leads to the consideration of individual integrals given by (2.5.1). These were then used in the context of 𝑊𝑜 given by a Mooney-Rivlin form. The M-R form is of special interest for sphere problems, since it can give either monotone or nonmonotone inflation response as determined by the M-R parameters and the shell thickness. Thus for example it was possible to show how residual stress could potentially cause a nonmonotone response to become monotone, or vice versa, without changing the problem geometry or the constitutive law. Here it is to be noted that the same integrals (2.5.1) will arise in an alternative treatment that retains 𝑊 = 𝑊𝑜 + 𝑊𝜏 with 𝑊𝜏 given by the first of (2.1.13) but with an alternative choice for 𝑊𝑜. More generally, it is immediate that the form 𝑊 = 𝑊𝑜+𝑊𝜏 with 𝑊𝜏 given by either of (2.1.13) has the property that the material response reduces to that of a conventional incompressible hyperelastic 49 solid with stored energy 𝑊𝑜 in the event that the residual stress 𝝉 vanishes. It is important to realize that a separate issue, and one that is not addressed in this paper, is whether or not a given non- vanishing 𝝉 could be associated with a specific relaxed configuration that is different from the reference configuration (meaning that there is no residual stress in the relaxed configuration). In other words, whether the hyperelastic response from the residually stressed reference configuration using 𝑊 = 𝑊𝑜 + 𝑊𝜏 matches the response from the relaxed configuration using 𝑊 = 𝑊𝑜 or perhaps a suitably modified 𝑊𝑜 (meaning a 𝑊𝑜 that does not incorporate a residual stress). Because this issue is not addressed here, care must be taken in describing the meaning of findings encapsulated in, for example, Figures 2.8 or 2.9. That is why our monotone vs. nonmonotone statement in the preceding paragraph avoids describing the results in terms of a residually stressed M-R material. Addressing these types of subtle issues gives rise to broader notions of how to frame such considerations in a fundamental way [54]. Within the specific context of Mooney-Rivlin type response these issues are addressed in [3], where there is a special focus on conditions of plane strain. Our attention here was also limited to residual stress fields with radial normal stress in the three-term form 𝜏𝑅𝑅 = 𝑘0 + 𝑘1𝑅𝑞1 + 𝑘2𝑅𝑞2 for suitably chosen integer exponents 𝑞2 > 𝑞1 > 0; these being the simplest nontrivial forms consistent with (2.2.9). However, such 𝜏𝑅𝑅 fields are of one sign, and the resulting 𝜏ΘΘ fields have exactly one internal node. For the purpose of constructing more general residual stress fields that are amenable to the treatment presented here, one may consider the general form containing 𝑛 terms (𝑛 ≥ 3): 𝜏𝑅𝑅 (𝑅) = 𝛼3𝑛−4 (cid:16) 𝑐1𝑅0 + 𝑐2𝑅2 + 𝑐3𝑅5 + ... + 𝑐𝑛−1𝑅3𝑛−7 + 𝑅3𝑛−4(cid:17) . (2.7.1) This form, with exponents given in the above particular way, is potentially highly workable for additional analysis because it generates integrals 𝑃 ¯𝜏|𝑞 given by (2.5.1) in terms of specific quotient expressions that do not make use of special functions. Note that the 𝜏𝑟𝑟 fields given by (2.3.3), (2.3.4) and (2.3.5) are all special cases of the form (2.7.1), but (2.3.1) is not. For example, (2.3.5) is 𝑛 = 5 with 𝑐2 = 𝑐3 = 0 and 𝑐1 and 𝑐4 given by complicated expressions. In general for (2.7.1), the parameter 𝛼3𝑛−4 can be regarded as an overall strength-of-field parameter, whereupon the degrees 50 of freedom conferred by the parameters 𝑐1 to 𝑐𝑛−1 is reduced by two so as to meet the boundary conditions (2.2.9). Thus the general form (2.7.1) effectively provides 𝑛-3 degrees of freedom beyond that of the overall strength-of-field for tuning the form of the residual stress fields. 51 CHAPTER 3 RESIDUAL STRESS FIELDS IN THE UNIFORM GROWING OF A HYPERELASTIC SPHERICAL SHELL Chapter 2 showed how the mechanical effect of the growth process in any later loading could be formally accounted for by making appropriate use of the associated residual stress field 𝝉 in the loading analysis. In particular, the examination in Chapter 2 revealed the sensitive dependence of subsequent mechanical loading on the specific features of 𝝉. With this insight, it is now useful to further inquire into the nature of the residual stress field 𝝉. To this end, following the same framework presented in [50], the treatment of growth-induced residual stress fields in a single-layer hyperelastic spherical shell is presented. Subsequently, the treatment is extended for a bilayer spherical shell utilizing this framework. One specific feature that requires clarification is the effect of the interface between gray and white matter on the residual stress field 𝝉. This requires a return to the broader kinematical considerations of the theory of finite growth elaborated in Section (1.3). The present chapter is organized as follows. The general hyperelastic treatment of a single layer spherical shell is given in Section (3.1), followed by the demonstration of the residual stress fields in normal directions for the considered five different growth cases in Section (3.2). The kinematic of volumetric growth of a bilayer spherical shell are introduced in Section (3.3). Key findings of this section, specifically the effect of mismatch growth of the layers on the residual stress fields, are presented in Section (3.4). In section (3.5), the findings of both analytical model is compared with the experimental studies introduced in Section (1.2.3). 3.1 Kinematics for volumetric growth of a single layer spherical shell We examine how the volumetric growth leads to residual stress fields within a spherical shell experiencing an inhomogeneous growth using the framework described in the previous section. To this end, a finite thickness spherical shell with inner radius 𝑅𝑖 > 0 and outer radius 𝑅𝑜 > 𝑅𝑖 prior to any loading or any growth is considered. The loading is taken to consist of applied pressures 𝑃𝑖 and 𝑃𝑜 on the inner and outer boundaries. These symmetric conditions motivate the consideration 52 of the radial deformation for volumetric growth 𝑟 = 𝑟 (𝑅), 𝜃 = 𝛩, 𝜙 = 𝛷 (3.1.1) on 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑜, 0 ≤ 𝛩 < 2𝜋, 0 ≤ 𝛷 ≤ 𝜋 where the radial inflation function 𝑟 (𝑅) is to be determined. Thus (3.1.1) is a map from reference spherical coordinates (𝑅, 𝛩, Φ) to deformed grown spherical coordinates (𝑟, 𝜃, 𝜙). Let {e𝑅, eΘ, eΦ} and {e𝑟, e𝜃, e𝜙} represent unit basis vectors in the spherical coordinate system of the respective reference, and deformed grown configurations. It follows from (3.1.1) that the total deformation gradient is given by F = 𝑟′(e𝑟 ⊗ e𝑅) + 𝜆 (cid:0)e𝜃 ⊗ eΘ + e𝜙 ⊗ eΦ (cid:1) , (3.1.2) with 𝑟′ = 𝑑𝑟/𝑑𝑅 and 𝜆 = 𝑟/𝑅. Here 𝜆 = 𝜆(𝑅) is the azimuthal stretch, meaning that it is the stretch along spherical surfaces. Because of the spherical symmetry, the stretch 𝜆 is the same in all directions upon each spherical surface. However 𝜆 may vary through the thickness of the shell. As detailed in previous section and specifically represented in (1.3.2), the kinematics of growth is expressed by two subsequent process: growth and elastic deformation. A spherically symmetric and continuous growth is considered. The growth tensor is given by Fg = 𝛾𝑟E𝑅 ⊗ E𝑅 + 𝛾𝜃 (EΘ ⊗ EΘ + EΦ ⊗ EΦ), (3.1.3) where {E𝑅, EΘ, EΦ} are the unit basis vectors for the reference B0 and virtual intermediate config- uration B𝜉. Here, 𝛾𝑟 = 𝛾𝑟 (𝑅) and 𝛾𝜃 = 𝛾𝜃 (𝑅) are the strictly positive growth functions in the radial and azimuthal directions respectively. Each growth function must be greater than one for actual volumetric growth in the given direction. Otherwise, it represents atrophy or resorption. Following the terminology used in [50], 𝛾𝑟 and 𝛾𝜃 are referred to as the radial and area growth respectively. Considering a spherical body, the radial growth increases the length of the shell in radial direction while the area growth increases the surface area. The case of 𝛾𝑟 = 𝛾𝜃 is called the isotropic growth, which is a straight forward dilatation of the body. The local volume change ratio associated with the growth following (1.3.4) is given by 𝑣 = 𝛾𝑟 𝛾2 𝜃 = 𝑣(𝑅). (3.1.4) 53 The functions 𝛾𝑟 (𝑅) and 𝛾𝜃 (𝑅) need not be continuous. If they are discontinuous then growth alone could generate discontinuities in the body. The elastic deformation tensor, which in the event of discontinuities growth serves to restore the continuity of the body, is described by Fe = 𝛼𝑟e𝑟 ⊗ E𝑅 + 𝛼𝜃 (cid:0)e𝜃 ⊗ EΘ + e𝜙 ⊗ EΦ (cid:1) , (3.1.5) in the deformed configuration B defined by the unit basis vector {e𝑟, e𝜃, e𝜙}, as illustrated in figure 1.4. Here, 𝛼𝑟 and 𝛼𝜃 are the elastic principal stretches in the radial and azimuthal directions respectively. The elastic incompressibility constraint (1.3.3) motivates the notational consolidation 𝛼 ≡ 𝛼𝜃, and 𝛼𝑟 = 𝛼−2. Thus, the form of (3.1.5) is revised, using (3.1.6), as Fe = 𝛼−2e𝑟 ⊗ E𝑅 + 𝛼(e𝜃 ⊗ EΘ + e𝜙 ⊗ EΦ), and the principal invariants are given, using (1.3.5), (1.3.9) and (3.1.7), by 𝐼1 = 𝛼−4 + 2𝛼2, 𝐼2 = 2𝛼−2 + 𝛼4. (3.1.6) (3.1.7) (3.1.8) For the deformation considered here, the combination of (1.3.2), (3.1.2), (3.1.3), and (3.1.7) gives that 𝑟 ′ = 𝛼−2𝛾𝑟, and 𝜆 = 𝛼𝛾𝜃. The volume constraint (1.3.4) becomes det F = 𝜆2 𝑟′ = 𝑟 2 𝑟′/𝑅2, making 𝑟 ′ = 𝑣 𝜆2 Also, the total deformation in the form 𝑟 2 𝑑𝑟 = 𝑣𝑅2 𝑑𝑅 integrates to 𝑟 3 = 𝑟 3 𝑖 + 3 ∫ 𝑅 𝑅𝑖 𝑣(𝜁) 𝜁 2 𝑑𝜁, (3.1.9) (3.1.10) (3.1.11) where 𝑟𝑖 = 𝑟 (𝑅𝑖) and 𝜁 is a dummy integration variable. More generally (3.1.11) provides the map 𝑟 = 𝑟 (𝑅) in terms of the single parameter 𝑟𝑖 which still needs to be determined. The Cauchy stress tensor takes the form 𝝈 = 𝜎𝑟𝑟 (e𝑟 ⊗ e𝑟) + 𝜎𝜃𝜃 (e𝜃 ⊗ e𝜃 + e𝜙 ⊗ e𝜙), (3.1.12) 54 with the combination of (1.3.10), (3.1.7) and (3.1.8) yields to 𝜎𝑟𝑟 = −𝑝 + 2𝛼−4 𝜕𝑊 𝜕𝐼1 𝜎𝜃𝜃 = −𝑝 + 2𝛼2 (cid:16) + 2 𝜕𝑊 𝜕𝐼1 , + 4𝛼−2 𝜕𝑊 𝜕𝐼2 𝛼−2 + 𝛼4(cid:17) 𝜕𝑊 𝜕𝐼2 . The equilibrium equation (1.3.7) gives that 𝑝 = 𝑝(𝑅) along with the requirement 𝑑𝜎𝑟𝑟 𝑑𝑟 + 2 𝑟 (𝜎𝑟𝑟 − 𝜎𝜃𝜃) = 0. (3.1.13) (3.1.14) (3.1.15) The specified pressures 𝑃𝑖 and 𝑃𝑜 at the inner and outer surfaces yield the boundary conditions 𝜎𝑟𝑟 (cid:12) (cid:12)𝑟𝑖 = −𝑃𝑖, 𝜎𝑟𝑟 (cid:12) (cid:12)𝑟𝑜 = −𝑃𝑜, (3.1.16) where 𝑟𝑜 = 𝑟 (𝑅𝑜) denotes the deformed outer radius of the spherical shell. One is easily led to the analytical expression to determine 𝑟𝑖 by performing the relevant substitutions of (3.1.13), (3.1.14), (3.1.10) into the equilibrium equation (3.1.15), integrating, and applying the boundary conditions (3.1.16). Once 𝑟𝑖 is known, the deformation 𝑟 = 𝑟 (𝑅) through (3.1.11) and the elastic deformation tensor Fe is fully specified. Such a treatment could be applied however the direct use of the stored energy density offers more straightforward procedure to obtain the output response of the shell to the growth. For the purpose of this study, a useful and elegant procedure is obtained by rewriting 𝑊 (𝐼1, 𝐼2) as 𝑤(𝛼) upon making use of (3.1.8). Then note that 𝜕𝑤 𝜕𝛼 = 𝜕𝑊 𝜕𝐼1 (cid:16) 𝜕𝐼1 𝜕𝑊 𝜕𝐼2 + 𝜕𝛼 𝜕𝐼2 𝜕𝛼 𝛼 − 𝛼−5(cid:17) 𝜕𝑊 𝜕𝐼1 (𝜎𝜃𝜃 − 𝜎𝑟𝑟) = = 4 = 2 𝛼 𝑟 𝛼 𝑑𝜎𝑟𝑟 𝑑𝑟 + 4 (cid:16) 𝛼3 − 𝛼−3(cid:17) 𝜕𝑊 𝜕𝐼2 , (3.1.17) where the last two steps employed first the pair (3.1.13) & (3.1.14) and then (3.1.15). The chain of equations in (3.1.17) causes the equilibrium equation in (3.1.15) to take the form 𝑑𝜎𝑟𝑟 𝑑𝑟 𝛼 𝑟 𝜕𝑤 𝜕𝛼 = 𝑑𝑟. 55 (3.1.18) At this point, with the aid of (3.1.9) and (3.1.10), developing the following derivation will be useful to derive an alternative form of (3.1.18), namely 𝑑𝛼 𝑑𝑟 = 𝑑𝛼/𝑑𝑅 𝑑𝑟/𝑑𝑅 = 𝑣 − 𝜆3 − 𝜆2𝛼𝑅𝛾′ 𝜃 𝑣𝑅𝛾𝜃 , (3.1.19) where 𝛾′ 𝜃 = 𝑑𝛾𝜃/𝑑𝑅. With the aid of the boundary conditions (3.1.16), and using (3.1.19), one transforms (3.1.18) into Δ𝑃 ≡ 𝑃𝑖 − 𝑃𝑜 = ∫ 𝛼𝑜 𝛼𝑖 𝑣 𝑣 − 𝜆3 − 𝜆2𝛼𝑅𝛾′ 𝜃 𝜕𝑤 𝜕𝛼 𝑑𝛼, (3.1.20) where 𝛼𝑖 = 𝛼(𝑅𝑖) = 𝑟𝑖/𝛾𝜃 (𝑅𝑖)𝑅𝑖 and 𝛼𝑜 = 𝛼(𝑅𝑜) = 𝑟𝑜/𝛾𝜃 (𝑅𝑜)𝑅𝑜 are the elastic stretches at the inner and outer surface respectively. Also note that 𝛼𝑜 can be expressed in terms of 𝛼𝑖. All together this shows that Δ𝑃 is indeed the function of 𝛼𝑖. Namely, in view of (3.1.9) and (3.1.11), the elastic stretches 𝛼𝑖 and 𝛼𝑜 are related by 𝛼3 𝑜 = (cid:18) 𝛾𝑟 𝛾𝜃 1 − (cid:19) 𝑅3 𝑖 𝑅3 𝑜 + 𝛼3 𝑖 (cid:19) 3 . (cid:18) 𝑅𝑖 𝑅𝑜 (3.1.21) Dividing the numerator and denominator of the above expression by 𝛾3 𝜃 it is convenient to define Γ = 𝛾𝑟/𝛾𝜃 = Γ(𝑅), and Ξ = −𝛼3𝑅𝛾′ 𝜃/𝛾𝜃 = Ξ(𝑅) whereupon (3.1.20), utilizing (3.1.4) and (3.1.9), takes the alternative form Δ𝑃 = ∫ 𝛼𝑜 𝛼𝑖 Γ Γ + Ξ − 𝛼3 𝜕𝑤 𝜕𝛼 𝑑𝛼, (3.1.22) which presents a general relation between the amount of growth, applied pressure, and deformation resulting from growth. Equation (3.1.22) gives the residual stress fields because of the growth in the absence of the loading. It also shows that the growth rate might affect the residual stress fields due to the existence of 𝛾′ 𝜃 in Ξ in addition to the growth functions [30]. In evaluating (3.1.22) both Γ and Ξ may vary with 𝑅 so the integration with respect to 𝛼 requires a parameterization of Γ and Ξ as a function of 𝛼 to properly account for 𝑣 within the integral, i.e., Γ(𝑅) = ˜Γ(𝛼), and Ξ(𝑅) = ˜Ξ(𝛼). 56 For the remainder of this Section 3.1, attention is confined to a special case: homogeneous growth of the spherical shell meaning that the growth functions 𝛾𝑟 and 𝛾𝜃 are independent of 𝑅. Thus Γ is a constant scalar, and Ξ = 0. Equation (3.1.22) yields to Δ𝑃 = ∫ 𝛼𝑜 𝛼𝑖 Γ Γ − 𝛼3 𝜕𝑤 𝜕𝛼 𝑑𝛼. (3.1.23) The classical incompressible isotropic hyperelastic result with no growth is achieved upon setting Γ = 1 (see for example (7.18) of [55] and (5.3.21) of [91]). This setting also corresponds to the isotropic growth (pure dilatation) that produces no residual stress in the shell. Equation (3.1.23) also retrieves the generalized swelling form that given in (28) of [150] and (22) of [97]. The form of constitutive model requires to be specified at this point so as to determine the stress fields. In this regard, we shall consider the well-known Mooney-Rivlin model to be used in specific examples, namely 𝑊MR(𝐼1, 𝐼2) = 1 2 𝜇 (cid:0)𝜅(𝐼1 − 3) + (1 − 𝜅) (𝐼2 − 3)(cid:1), (3.1.24) where 𝜇 is the shear modulus with the unit of stress. 𝐼1 and 𝐼2 are the principal invariants given already in (3.1.8). 𝜅 is a scalar material parameter in the range 0 ≤ 𝜅 ≤ 1 to be consistent with the Baker-Ericksen inequalities. The shear modulus 𝜇 vary with the position within the spherical shell, i.e., 𝜇 = 𝜇(X); however, we restrict attention to homogeneous shear modulus for now despite of the possibility of growth-induced changes. One enables to obtain the auxiliary form of the strain energy function, using (3.1.8), that 𝑤MR(𝛼) = 1 2 𝜇 (cid:0)𝜅(𝛼−4 + 2𝛼2 − 3) + (1 − 𝜅) (2𝛼−2 + 𝛼4 − 3)(cid:1). (3.1.25) 3.2 Growth-Induced residual stress fields in single layer spherical shell For the case of material model (3.1.24), the combination of (3.1.23) and (3.1.25) yields to an integral that provides a relation between: amount of growth, applied pressure Δ𝑃, and the deformation of inner surface. By performing integration, and substitution of the boundaries in the view of (3.1.21), one can readily obtain an analytical expression where the parameter 𝑟𝑖 is the only unknown. In the absence of pressurization (Δ𝑃 = 0), the stress field corresponds to the 57 Figure 3.1 Radial and circumferential residual stresses in a growing shell 𝜎𝑟𝑟 (red), and 𝜎𝜃𝜃 (blue) generated by radial growth (left panel) and area growth (right panel). Here, 𝜇 = 1, 𝑅𝑖 = 1, 𝑅𝑜=2 with the following growth constants: 𝛾𝑟 = 2 and 𝛾𝜃 = 1 in the left panel and 𝛾𝑟 = 1 and 𝛾𝜃 = 2 in the right panel. √ residual stress field as a result of the growth. Thus the deformed inner radius 𝑟𝑖 can be determined by applying the boundary conditions and subsequently the whole kinematic becomes fully known through (3.1.11). Eventually, the radial residual stress 𝜎𝑟𝑟 and the circumferential stress 𝜎𝜃𝜃 can be determined via (3.1.22) and the equilibrium equation (3.1.15) respectively. Suppose now that the amount of growth is chosen in the way that the volume of the shell becomes double, i.e., 𝑣 = 2. Here we consider two alternative growth forms to make a quantitative characterization of the residual stress fields. Setting 𝛾𝑟 = 1 and 𝛾𝜃 = 1, equation (3.1.4) leads to √ 𝛾𝜃 = 2 and 𝛾𝑟 = 2. The growth constant pairs correspond to area (𝛾𝑟 = 2, 𝛾𝜃 = 1) and radial (𝛾𝑟 = 1, 𝛾𝜃 = √ 2) growth in the shell respectively. The left and right panel of figure 3.1 shows the radial and circumferential residual stress fields on account of radial and area growth respectively. Figure 4 also exhibits the residual stress fields for the Mooney-Rivlin form with alternative values of 𝜅. Figure 3.1 confirms that the growth causes the residual stresses even in the absence of any external loading. The same stress fields have already shown in Fig. 15.5 of [50] but for only 𝜅 = 1. It also shows that the different growth forms, i.e., radial and area growth, influence the qualitative behavior of the residual stress fields. Namely, the radial growth leads to compressive stress whereas 58 the area growth causes tensile stress along the radius of the spherical shell. The radial stress 𝜎𝑟𝑟 is of a single sign in both cases. However, the sign of 𝜎𝜃𝜃 coincides with the sign of 𝜎𝑟𝑟 in the inner surface, then the sign changes as goes to the outer surface of the shell. Hence, the circumferential stress 𝜎𝜃𝜃 can become either compressive or tensile depending on the location. The residual stress profiles produced by area and radial growth have also different characteristics. Most notable is that the radial growth shows a localization tendency towards the inner surface whereas the area growth leads to a more evenly distributed stress response. The value of 𝜅 slightly affect the qualitative behavior of the residual stress regardless of the form of growth. Although the qualitative behavior stays the same for all cases, the magnitude of residual stress increases as the value of 𝜅 decreases. Our focus is here to understand the behavior of residual stress under different growth conditions instead of a particular response of a material model. Therefore, the attention is restricted for developing a treatment for the following examples using 𝜅 = 1 that correspond to the well-known Neo-Hookean constitutive model 𝑊NH(𝐼1) = 𝜇 2 (𝐼1 − 3). (3.2.1) Considering the spherical symmetry, (3.2.1) can be written in the form of ˜𝑊 (𝛼𝑟, 𝛼𝜃) = 1 2 𝜇(𝛼2 𝑟 + 2𝛼2 𝜃 − 3) then one enables to obtain, using (3.1.6), that 𝑤(𝛼) = 𝜇 2 (𝛼−4 + 2𝛼2 − 3). For the case of material model (3.2.2), (3.1.23) turns in the following form Δ𝑃 = 2𝜇Γ ∫ 𝛼𝑜 𝛼𝑖 𝛼 − 𝛼−5 Γ − 𝛼3 𝑑𝛼. (3.2.2) (3.2.3) For the purposes in what follows, the section continues to provide the residual stress fields for alternative growth cases by applying the same solution procedure to (3.2.3). The alternative growth forms are determined as follows. Following the aforementioned specific example, the growth field is again set as a constant and for 𝑣 = 2. The radial and area growth forms showed in the previous example are included and labeled as the case I and V in table 3.1. Three alternative pairs of growth constants are considered along with these two forms. As the third case, the isotropic growth where 59 𝛾𝑟 = 𝛾𝜃 thus Γ = 1 is considered. Eventually, two more cases for the intermediate Γ values are generated by taking into consideration the growth field constraint. Table 3.1 Five alternative growth constants pairs all giving 𝑣 = 2. Case 1 2 3 4 5 𝛾𝑟 2 1.65 √ 3 2 1.13 1 𝛾𝜃 1 1.10 √ 3 2 1.32 √ 2 Γ 2 1.5 1 0.85 0.70 Type of Growth radial growth only radial growth > area growth isotropic growth area growth > radial growth area growth only Figure 3.2 Residual stress fields with 𝜎𝑅𝑅 in the left panel, and 𝜎𝜃𝜃 in the right panel. Here, 𝜇 = 1, 𝑅𝑖=1, 𝑅𝑜=2 and with the following growth constants ratio given in the table 3.1 as Γ = 2 (solid), Γ = 1.50 (solid with square marker), Γ = 1 (dotted), Γ = 0.83 (dash-dotted), and Γ = 0.707 (dashed). Plots of 𝜎𝑟𝑟 and 𝜎𝜃𝜃 shown in figure 3.2 reflects how the type of growth affect the behavior of residual stress fields. The sign of residual stress is consistent with the observation we made in the previous example. The sign of 𝜎𝑟𝑟 depends on whether 𝛾 is less than or greater than one. Namely, the cases in which radial growth is dominant or Γ > 1 generates a compressive residual stress field whereas the cases in which Γ < 1 leads to the tensile compressive stress in the radial direction. The radial stress vanished at the boundaries and is of single sign in for all cases. In contrast, the circumferential residual stress takes both positive and negative values for all cases. The sign of 𝜎𝜃𝜃 is initially the same with the sign of 𝜎𝑟𝑟 associated with the same growth type 60 but then changes towards the outer radius. The location where the sign change closer to the outer radius as Γ decreases. 𝜎𝑟𝑟 show a nonmonotonic behavior involving two intervals: increase to a local maximum followed by a monotone decrease with 𝑅 → 𝑅𝑜. However, 𝜎𝜃𝜃 is one of continued monotone the magnitude of residual stress decreases as 𝑅 → 𝑅𝑜. Therefore, the maximum residual in radial direction stress occurs somewhere closer to the inner surface whereas it occurs at the inner boundary in the circumferential direction. Differently, if the growth is isotropic (Γ = 1), the residual stress field in both radial and circumferential direction is vanished for everywhere. Figure 3.3 The influence of growth on the inflation response as determined by (3.2.3) for five different growth constant pairs given in table 3.1, 𝑅𝑖=1, 𝑅𝑜=2 and 𝜇=1. Lastly, we shall provide a demonstration of the pressure-inflation response of a growing spherical shell for the growth constants pairs given in table 3.1. One can also determine the pressure-inflation response by applying the aforementioned treatment to (3.2.3). This gives plots of Δ𝑃 with respect to the elastic stretch at the inner boundary 𝛼𝑖, as shown in figure 3.3. The negative portion 61 of pressure-inflation graphs is not displayed. Following the standard classification scheme for spherical inflation [27], all cases show a so-called type (b) behavior, meaning that Δ𝑃 increase to a local maximum followed by monotone decrease. Positive pressurization (Δ𝑃 > 0) leads to inflation for all cases with the different local maximum values ranging between 1.73 and 0.56. In the absence of the pressurization Δ𝑃 = 0, isotropic growth leads to dilatation (𝑟 3 = 𝑣𝑅3) that gives 𝑟 ′ = 𝜆. Here, 𝛼𝑎 is varying from 0.74 (Case-I) to 1.15 (Case-V) depending on the type of growth. It can be conclude that as the area growth gets dominated the inner radius of the sphere is pushed out more, depending on the Γ. In other words, smaller Γ values correspond to higher strain levels at the inner surface of the sphere. 3.3 Kinematics for volumetric growth of a bilayer spherical shell In Section 3.1, the residual stress fields in the single-layer Neo-Hookean type spherical shell for the alternative growth forms are presented. The treatment of single-layer shells provided useful insight into the residual stress fields associated with the growth field. However, the multi-layered structures are more predominant among the biological tissues and organs, for instance, arteries [56], airways [52], and brain [106]. In case of multilayered structures, the relative stiffness, growth rate and thickness ratio of each layer could play a role in the behavior of residual stress fields. One of the hypothesis, which explains the cortical folding, already based on the different tangential growth ratio of gray matter (top layer) and white matter (inner layer). A recent experimental study, which presents the residual stress field in a developing ferret brain, also supports this hypothesis [144]. The differential growth conditions need to be addressed to understand the residuals stresses in a growing multi-layered structure. A treatment is therefore developed to determine the qualitative behavior of residual stress fields in a bilayer spherical shell following the framework used in section 3.1. The reference configuration is now partitioned into two parts: an inner core 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚 and an outer shell 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜 where 𝑅𝑚 denotes the initial interface radius. Each experiences homogeneous growth but with possibly different values for 𝛾𝑟 and 𝛾𝜃 62 𝛾𝑟 =    𝛾𝑖𝑟 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚 𝛾𝑜𝑟 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜 𝛾𝜃 =    𝛾𝑖𝜃 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚 𝛾𝑜𝜃 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜 (3.3.1) where the 𝑖 and 𝑜 subscripts stands for inner and outer layers respectively. Similar to (3.1.4), the local volume change ratio of layers can be expressed as 𝑖𝜃, 𝑣𝑖 = 𝛾𝑖𝑟 𝛾2 and 𝑣𝑜 = 𝛾𝑜𝑟 𝛾2 𝑜𝜃. (3.3.2) where 𝑣𝑖 and 𝑣𝑜 represent the total change in volume of each layer. Note that 𝑣𝑖 and 𝑣𝑜 are constant as a function of 𝑅 due to homogeneous growth. The relations given in (3.1.9) and (3.1.10) can be applied to the inner and outer layer individually with an easy subscript modification. It is assumed that the inner and outer are perfectly bonded. Integrating (3.1.10) from the interface to a generic radial value 𝑅 gives ∫ 𝑅𝑚 𝑟 3(𝑅) = 𝑟 3 𝑚 − 3 𝑟 3 𝑚 + 3    𝑅 ∫ 𝑅 𝑅𝑚 𝑣𝑖 (𝜁) 𝜁 2𝑑𝜁, 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚, 𝑣𝑜 (𝜁) 𝜁 2𝑑𝜁, 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜, (3.3.3) where 𝑟𝑚 = 𝑟 (𝑅𝑚). Due to the compatibility condition, there is no need to use the separate notation for the parameter 𝑟𝑚. The deformation mapping (3.3.3) enables one to obtain the following form of the radial deformation for each layer 𝑟 (𝑅) =    (cid:0)𝑟 3 𝑚 − 𝑣𝑖 (𝑅3 𝑚 − 𝑅3)(cid:1) 1/3, 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚 (cid:0)𝑣𝑜 (𝑅3 − 𝑅3 𝑚) + 𝑟 3 𝑚 (cid:1) 1/3, 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜. (3.3.4) Here (3.3.4) provides a map for the deformation in terms of 𝑟𝑚 which is the only unknown parameter. With the aid of the volume constraint (3.1.9), the elastic stretches at the inner and outer layer are written, using (3.3.4), as 63 𝛼3(𝑅) = (cid:19) 3 (cid:18) 𝑟𝑚 𝛾𝑖𝜃 𝑅 − Γ𝑖 (cid:18) 𝑅3 𝑚 𝑅3 (cid:19) − 1 , 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚, (cid:19) 3 (cid:18) 𝑟𝑚 𝛾𝑜𝜃 𝑅 + Γ𝑜 (cid:18) 1 − (cid:19) 𝑅3 𝑚 𝑅3 , 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜, (3.3.5)    where Γ𝑖 = 𝛾𝑖𝑟/𝛾𝑖𝜃 and Γ𝑜 = 𝛾𝑜𝑟/𝛾𝑜𝜃. Following the derivations from (3.1.15) to (3.1.23) in the previous section, one can derive the general form of the stress equation for a spherical deformation in which a homogeneous growth is taken into account. Similarly, in case of a bilayer spherical shell, (3.1.23) generalizes to Δ𝑃 ≡ 𝑃𝑖 − 𝑃𝑜 = ∫ 𝛼𝑟𝑚 𝛼𝑟𝑖 Γ𝑖 Γ𝑖 − 𝛼3 𝜕𝑤(𝛼) 𝜕𝛼 𝑑𝛼 + ∫ 𝛼𝑟𝑜 𝛼𝑟𝑚 Γ𝑜 Γ𝑜 − 𝛼3 𝜕𝑤(𝛼) 𝜕𝛼 𝑑𝛼. (3.3.6) where 𝛼𝑟𝑖 = 𝛼(𝑅𝑖), 𝛼𝑟𝑜 = 𝛼(𝑅𝑜) and 𝛼𝑟𝑚 = 𝛼(𝑅𝑚) are the elastic stretches, already given in (3.3.5), at the inner, outer and interface respectively, and Δ𝑃 denotes the pressure difference across the complete two layer spherical shell. Equation (3.3.6) provides a relation between the amount of growth of each layer, the applied pressure and the deformed interface radius 𝑟𝑚. In the absence of pressurization i.e., Δ𝑃 = 0, (3.3.6) corresponds the residual stress fields in the layers due to growth. At this point, the material model should be specified to obtain the qualitative behavior of the residual stress field. A neo-Hookean material model given in (3.2.1) but with a possibly different shear modulus 𝜇 is utilized for each individual layer 𝜇 =    𝜇𝑖, 𝑅𝑖 ≤ 𝑅 ≤ 𝑅𝑚, 𝜇𝑜, 𝑅𝑚 ≤ 𝑅 ≤ 𝑅𝑜. (3.3.7) For the case of material model (3.2.1) with layer-specific shear moduli (3.3.7), in the absence of pressurization (3.3.6) takes to form of 𝜇𝑖Γ𝑖 ∫ 𝛼𝑟𝑚 𝛼𝑟𝑖 𝛼 − 𝛼−5 Γ𝑖 − 𝛼3 𝑑𝛼 + 𝜇𝑜Γ𝑜 ∫ 𝛼𝑟𝑜 𝛼𝑟𝑚 𝛼 − 𝛼−5 Γ𝑜 − 𝛼3 𝑑𝛼 = 0, (3.3.8) in which enable one to obtain the relation between the amount of growth and the deformation of neo- Hookean spherical shell. Equation (3.3.8) can be seen the additive form of the separate responses of the inner and outer layer. The explicit integration of (3.3.8) gives an analytical expression of 64 a function of 𝑟𝑚 which is the single unknown parameter. Once 𝑟𝑚 is determined from solving (3.3.8), the deformation of the inner and outer layer is obtained through (3.3.4). Subsequently, the residual stress in radial direction 𝜎𝑟𝑟 can be determined using the integration associated with the layer. Eventually, 𝜎𝜃𝜃 can be determined from the equilibrium equation (3.1.14) for each layer. While (3.3.8) has here been obtained as a formal process based on (3.3.6) for the case Δ𝑃 = 0, it is also useful to note its physical meaning. The situation at hand accounts to the consideration of boundary conditions (3.1.16) with 𝑃𝑖 = 0 and 𝑃𝑜 = 0. The radial normal stress 𝜎𝑟𝑟 continues to obey (3.1.18) in each layer whereupon (3.1.18) may be integrated across each layer making use of the substitutions that led from (3.1.18) to (3.1.23). Let 𝑟 − 𝑚 = 𝑟 (𝑅− 𝑚) and 𝑟 + 𝑚 = 𝑟 (𝑅+ 𝑚) denote locations on the interface with 𝑟 − 𝑚 on the core side of the interface and 𝑟 + 𝑚 on the shell side. It then follows that and 𝜎𝑟𝑟 (𝑟 − 𝑚) = ∫ 𝛼𝑟𝑚 𝛼𝑟𝑖 Γ𝑖 Γ𝑖 − 𝛼3 𝜕𝑤 𝜕𝛼 𝑑𝛼, 𝜎𝑟𝑟 (𝑟 + 𝑚) = ∫ 𝛼𝑟𝑜 𝛼𝑟𝑚 Γ𝑜 Γ𝑜 − 𝛼3 𝜕𝑤 𝜕𝛼 𝑑𝛼, whereupon (4.8) is simply the statement that 𝜎𝑟𝑟 (𝑟 − 𝑚) = 𝜎𝑟𝑟 (𝑟 + 𝑚) (3.3.9) (3.3.10) (3.3.11) for the neo-Hookean model. Because there are no shear tractions operating on 𝑟 = 𝑟𝑚 the formal analysis has recovered the required condition of traction continuity at the interface. The stress component 𝜎𝜃𝜃 does not enter into the determination of the fractions at 𝑟 = 𝑟𝑚 and so need not be continuous at the interface. Indeed it follows from (3.1.15) and (3.3.11) that 𝜎𝜃𝜃 (𝑟 + 𝑚) − 𝜎𝜃𝜃 (𝑟 − 𝑚) = 𝑟𝑚 2 (cid:18) 𝑑𝜎𝑟𝑟 𝑑𝑟 (cid:12) (cid:12) (cid:12)𝑟 + 𝑚 − 𝑑𝜎𝑟𝑟 𝑑𝑟 (cid:12) (cid:12) (cid:12)𝑟 − 𝑚 (cid:19) . (3.3.12) Consequently discontinuity in 𝜎𝜃𝜃 at 𝑟 = 𝑟𝑚 correlates with discontinuity in the slope of 𝜎𝑟𝑟 at 𝑟 = 𝑟𝑚. 65 3.4 Growth-Induced residual stress fields in bilayer spherical shell Figures 3.4, 3.5, 3.6, 3.7, and 3.8 display the residual stress field in the bilayer spherical shell for all possible combinations of alternative growth fields in the table 3.1. The left and right panels of these figures show the radial 𝜎𝑟𝑟 and circumferential residual stress 𝜎ΘΘ fields, respectively. The thickness of layers is considered equal to each other; the location of the interface 𝑅𝑚 = 1.5 is shown by a dashed black line for all cases. It is assumed that both layers have same shear modulus (i.e., 𝜇𝑖 = 𝜇𝑜 = 𝜇) in the cases presented in Figs (3.4 - 3.8) and both stress field 𝜎𝑟𝑟 and 𝜎𝜃𝜃 is normalized by 𝜇. The radial stress field 𝜎𝑟𝑟 is continuous for all cases and vanished only at the boundaries. Except for the pairs of growth fields (I-I), (I-II), (I-III), (I-IV), (II-I), (II-II) and (II-III), 𝜎𝑟𝑟 has a monotonically either increasing or decreasing behavior. The radial residual stress field tends to localize to the inner radius 𝑅 = 𝑅𝑖 as the radial growth is more dominant within the inner layer as shown in figure 3.4 and 3.5. In contrast, the circumferential normal stress 𝜎𝜃𝜃 is discontinuous except for the cases where both inner and outer layers undergo the same growth field. Regardless of the type of growth, the differential growth leads to a stress jump at the interface in the circumferential direction. Lastly, the behavior of circumferential stress fields remain the same in the outer layer although the inner layer undergoes different growth fields. It has monotonically increasing or decreasing behavior for all cases. The magnitude of stress field shifts in the vertical direction gradually as the growth field became area growth dominated Γ → 0.707. That results in a transition from the compressive to the tensile stress field in the tangential direction. 66 Figure 3.4 Residual stress fields 𝜎𝑟𝑟 and 𝜎𝜃𝜃 for a growing bilayer spherical shell. The growth field ratio of inner layer is Γ = 2 while the outer layer undergoes a growth with following ratio: Γ = 2 (solid), Γ = 1.5 (dashed), Γ = 1 (dash-dotted), Γ = 0.85 (dotted), and Γ = 0.70 (solid with circle marker). Figure 3.5 Residual stress fields 𝜎𝑟𝑟 and 𝜎𝜃𝜃 for a growing bilayer spherical shell. The growth field ratio of inner layer is Γ = 1.5 while the outer layer undergoes a growth with following ratio: Γ = 2 (solid), Γ = 1.5 (dashed), Γ = 1 (dash-dotted), Γ = 0.85 (dotted), and Γ = 0.70 (solid with circle marker). As displayed in previous figures clearly, the difference in growth fields of the layers leads to the stress discontinuity in the circumferential direction at the interface. Figure 12 exhibits the circumferential stress difference for each case presented above. The continuous cases in which the inner and outer layers have the same growth field are not displayed. The red squares show a transition from compressive to tensional while the blue circles represent a transition from tensional 67 Figure 3.6 Residual stress fields 𝜎𝑟𝑟 and 𝜎𝜃𝜃 for a growing bilayer spherical shell. The growth field ratio of inner layer is Γ = 1 while the outer layer undergoes a growth with following ratio: Γ = 2 (solid), Γ = 1.5 (dashed), Γ = 1 (dash-dotted), Γ = 0.85 (dotted), and Γ = 0.70 (solid with circle marker). Figure 3.7 Residual stress fields 𝜎𝑟𝑟 and 𝜎𝜃𝜃 for a growing bilayer spherical shell. The growth field ratio of inner layer is Γ = 0.85 while the outer layer undergoes a growth with following ratio: Γ = 2 (solid), Γ = 1.5 (dashed), Γ = 1 (dash-dotted), Γ = 0.85 (dotted), and Γ = 0.70 (solid with circle marker). to a compressive residual stress field. The black lines are obtained by linear fit in order to illustrate the change in the stress difference corresponding growth field of the outer layer. The highest slope is observed in the case where the inner layer grows purely in the radial direction. 68 Figure 3.8 Residual stress fields 𝜎𝑟𝑟 and 𝜎𝜃𝜃 for a growing bilayer spherical shell. The growth field ratio of inner layer is Γ = 0.70 while the outer layer undergoes a growth with following ratio: Γ = 2 (solid line), Γ = 1.5 (dashed line), Γ = 1 (dash-dotted), Γ = 0.85 (dotted), and Γ = 0.70 (solid with circle marker). The influence of relative shell thickness on the stress field is now investigated for only the case of the pair of (IV-II). Figure 3.10 shows the circumferential stress field 𝜎𝜃𝜃 along the normalized undeformed radius ¯𝑅 as the location of interface is given by the three alternatives 𝑅𝑚 = 1.25, Figure 3.9 Circumferential stress difference at the interface with respect to the growth ratio of the inner and outer layer. 69 Figure 3.10 The effect of thickness of inner and outer layer on the circumferential stress difference with respect to normalized radius for the case of the pair of (IV-I). Here 𝑅𝑖 = 1, 𝑅𝑜 = 2, 𝜇𝑖 = 𝜇𝑜 = 1 and with a) 𝑅𝑚 = 1.25 (dashed), b)𝑅𝑚 = 1.50 (solid), 𝑅𝑚 = 1.75 (dotted). 𝑅𝑚 = 1.50, and 𝑅𝑚 = 1.75 (all while keeping 𝑅𝑖 = 1 and 𝑅𝑜 = 2 ). The circumferential stress difference is illustrated by blue arrow for the case 𝑅𝑚 = 1.25 and denoted by Δ𝜎𝜃𝜃. As the inner layer thickness decreases, the behavior of the stress field of this layer transforms monotonically decreasing into monotonically increasing. On the other hand, only the slope of the stress field and the magnitude of stress at the interface change in the outer layer. For two subcases 𝑅𝑚 = 1.50 and 𝑅𝑚 = 1.75, the sign of stress field is different for either side of interface. However, the residual stress field has negative sign (compressive) at both side of the interface for the case of 𝑅𝑚 = 1.25. This sign transition is actually occur at 𝑅𝑚 = 1.37. The residual stress at the outer boundary 𝑅 = 𝑅𝑜 is the same for all cases. One also observes that the magnitude of stress difference at the interface gradually increases as the inner layer becomes thinner. The residual stress fields of the bilayer sphere with growing dependent material stiffness is illustrated by considering the same the growth field pair (IV-I) in figure 3.11. For this purpose, two alternative cases are considered: relatively stiff outer layer 𝜇𝑖/𝜇𝑜 = 1/10 and relatively stiff inner layer 𝜇𝑖/𝜇𝑜 = 10 in addition to the identically stiff both layers 𝜇𝑖/𝜇𝑜 = 1 that corresponds to blue 70 dashed curve in figure 3.7. The thickness ratio of layers is considered equal meaning that 𝑅𝑚 = 1.5 whose location is showed by dashed blue curve in both panels. Similarly, the radial stress field is continuous and vanished at the inner and outer boundaries. The behavior of radial stress field is the same which is monotonically decreasing regardless of the relative stiffness ratio. The relatively stiff inner layer leads to non-monocity specifically parabolic profile whereas a monotonically decreasing behavior is observed for the relatively compliant inner layer. Figure 3.11 The effect of relative shear modulus on the residual stress fields for the case of pair (IV-I) with three alternative shear modulus ratio: 𝜇𝑖/𝜇𝑜 = 1 (solid), 𝜇𝑖/𝜇𝑜 = 1/10 (dotted), 𝜇𝑖/𝜇𝑜 = 10 (dashed). 3.5 Comparison of results with experimental findings The number of experimental studies regarding the stress fields in the cortex is highly limited. However, there are some studies providing an understanding of the formation of growth-induced residual stress fields within the cortex. The objective of this section is to critically assess the results obtained from the single and bilayer hyperelastic spherical models in the context of the available experimental findings. As detailed in section 1.2.3, there are two experimental works presenting the residual stress fields in adult mouse brain [143] and evolving ferret brain [144]. The main distinction between these two works undeniably lies in the anatomical characteristics of the experimental subjects. The mouse brain is lissencephalic meaning that the outer surface of the brain is smooth. In contrast to 71 11.21.41.61.82-10-5051011.21.41.61.82-8-6-4-202 species with larger brain sizes such as humans, the mouse brain does not experience a gyrification process. Therefore, the experiment conducted on the mouse brain allows us to make a series of qualitative observations regarding the stress field in that cortical folding does not play any role. It is also known that the human brain experiences radial and circumferential growth until the 26th week of gestation [31]. Hence, the experimental findings could also be useful to understand the mechanical state of tissue just before the onset of cortical folding. In the mouse dissection experiment [144], only the radial cuts at various depths have been made to investigate the stress field in the coronal slice of the mouse brain. The cuts relieve the stress field perpendicular to the cutting direction. Hence, it can give insight about only the circumferential residual stress field 𝜎𝜃𝜃 during the analysis because the experimental data is available in this direction. The cutting locations and their depths have been shown in figure 1.2 by solid and open arrowheads. Figure 1.2 allows us to estimate the qualitative behavior of residual stress in the gray matter (GM) (Cf. figure 1.2 a2-b2), gray-white matter (GM-WM) interface (Cf. figure 1.2 a3-b3). The radial cut only through gray matter stayed closed, which means the compressional residual stress field across the cut. The relatively wider opening is observed at the GM-WM interfaces. The opening remains closed in the GM whereas the opening diminishes as goes toward the bulk of the tissue. It is therefore safe to assume that the white matter, starting from the interface, exhibits a tensional residual stress field characterized by a decreasing profile extending from the interface to the inner depths of the tissue. It is also noted that the openings become wider after 15 min the cuts were made as shown in figure 1.2 a4-b4, because of the viscoelastic nature of the brain soft tissue. Revisiting the results presented in section 3.2, figure 3.2 showed that both radial and area growth induce a residual stress field in circumferential direction, including both positive and negative values. Specifically, cases 4 and 5 given in Table 3.1, where the area growth is more dominant, lead to a transformation from compressive to tensile residual stress field as goes to outer radius to inner radius. This is similar to what was observed in the mouse brain experiment. It is worth noting that the analytical model predicts a stress field that exhibits continuity along the circumferential dimension. In the event that the stress profile remains continuous across the 72 interface, the circumferential stress profile should necessarily attain a zero value at some point due to the transition from a compression to a tensile stress field when traversing from the cortex to the subcortex. It is also reasonable to assert that the circumferential stress field should exhibit discontinuities at the interface. These discontinuities may be attributed to disparities in material properties or growth profiles between the outer and inner layers. In Section 3.4, a bilayer spherical model characterized by differential growth conditions was examined. This may more closely resemble the anatomical features of the mouse brain. The obtained results reveal that the stress field in the radial direction remains continuous, while the circumferential residual stress field exhibits discontinuities in scenarios where growth conditions differ between the layers. Recalling the stress field observed in the mouse brain, it is expected that the circumferential stress pattern 𝜎𝜃𝜃 should exhibit a positive sign in the region defined by 𝑅𝑖 < 𝑅 < 𝑅𝑚, followed by a transition to a negative sign in the region delimited by 𝑅𝑚 < 𝑅 < 𝑅𝑜. Consequently, it is reasonable to assert that the cases denoted as (IV-V), (III-V), (III-IV), (II-III), (II-IV), (II-V), (I-III), (I-IV), and (I-V) are likely to leads to a circumferential stress pattern akin to that observed in the mouse brain dissection experiment. Assuming both layers undergo an anisotropic growth process 𝛾𝑟 ≠ 𝛾𝜃, the resulting residual stress field 𝜎𝜃𝜃 emerges as a consequence of the outer layer experiencing more tangential growth. 73 CHAPTER 4 MICROMECHANIC INVESTIGATION OF THE GYRI-SULCI 4.1 Motivation and background Traumatic brain injury (TBI) is a broad term referring the pathologies in the brain that occur as a result of exogenous mechanical forces such as a blast or impact. The consequences of TBI, which is a major cause of long-term disabilities in both developed and developing countries, are dramatic - a recent study reported that an estimated ten million people will be affected annually by TBI [68]. Considering the Russian invasion of Ukraine, this number is likely to be higher than estimated. TBI leads to a wide range of impairments in physical and cognitive abilities that emerge either immediately (primary) or subsequently (secondary). Primary injuries (PI) are mostly associated with the transduced energy from the source of impact, which causes an abrupt injury, whereas the secondary injuries (SI) can be initiated via the PI or observed as a consequence of PI. The SI causes relatively long-lasting effects [25]. Along with the damage at the macro scale, physical forces lead to disruption of structural integrity, and dysfunction in the cellular and molecular processes within the neuron and glial cells. Due to the anisotropic and heterogeneous nature of the brain, the components of the central and peripheral nervous system can respond differently to the same intracranial stress and strain [61, 102], which causes distinctive impairments because of the specialization of the cells to perform a set of specific roles. An acute TBI can also evolve into progressive neuropathology such as Chronic Traumatic Encephalopathy (CTE). CTE is defined as "a distinct perivascular accumulation of hyperphospho- rylated tau (p-tau) in neurons and astrocytes within cerebral sulci" [113]. Tau is a protein that maintains the stabilization of the microtubule. The phosphorylated tau binds normal tau proteins and starts accumulating into neurofibrillary tangles (NFTs) [69]. The accumulation of p-tau in NFTs results in the change in morphology of neuron dendrite and neurotransmitter receptor expression, which is the main cause of impairments in cognitive abilities and memory loss. NFTs are the most common biomarkers of tauopathies such as CTE, and Alzheimer’s Disease (AD). The distribution of NFTs can diverge within the brain tissue depending on the neurodegenerative diseases. For 74 instance, NFTs are localized in the base of sulci and perivascular regions of the cortex in CTE [89]. The risk of developing CTE is higher in people exposed to the repetitive physical impact such as service members, athletes, hockey players, college, and professional football players. An MRI-based study showed that the repetitive force impact can lead to a change in the appearance of subtle hemorrhagic in the base of the sulcus, increases in the width of sulci in the frontal and occipital cortices, and decrease in total brain volume [78]. It should be noted that not only repetitive exposure, but a single moderate or severe TBI, could also suffice to initiate a tau pathology like CTE [132]. The initiation and progression of pathophysiological changes in a living subject is difficult to detect using currently available diagnostic tools [1]. Yet, histological data clearly show the distinctive p-tau pathology of CTE in the base of sulci and the spread of NFTs from the epicenter through white matter, as shown in figure 4.1. Figure 4.1 Distinct pattern of neurofibrillary tangles seen in CTE: at the depth of sulci (P-R and U-W), and around the small blood vessels (S, T, X, Y). (Reused from [84] page 47. In the public domain). Over the last decade the awareness regarding the long-term neurotraumatic consequences of TBI has been increased. To this end, several physical mechanisms have been proposed to explain mechanical deformation associated with the occurrence of CTE. One of the possible mechanisms is the strain and strain rate localizing at the depth of sulci. Computational models predicted that the depth of sulci is likely to experience a higher strain and strain rate during TBI [39], [47], [154], summarized in table 4.1. Regardless of the amplitude and duration of impact, the maximum strain and strain rate have been observed within the sulcus [47]. The author has also previously reported that the deviatoric stress is localized at the depth of the sulcus during blast-induced TBI [145]. The 75 other injury mechanism is proposed by Kornguth et al. [78], which is called "the water hammer effect". The water hammer effect is based on the idea that the exogenous impact is transmitted to the base of the sulci. CSF is driven into the sulci as the elastic brain moves towards the relatively rigid calvarium and non-compressible CSF. Thus, the base of the sulcus is exposed to the impulse being transmitted from the skull to CSF, which leads to the damage specifically at the interface of gray and white matter. Table 4.1 The details of the computational models simulating the gyri-sulci micromechanics. Single sulcus Single sulcus Study Domain [16] [76] [124] Brain (3D) [39] Injury Mechanism Rigid impactor Cavitation bubble expansion Blast exposure (Shock-Tube) Axial Brain Slice (2D) Rigid wall impact [47] Brain (3D) Helmet-to-helmet impact Fall impact Road accident Strain/Strain Rate 0.2-0.45/- 0.41/- 0.070-0.110 / 5-12 s−1 0.1-0.25 / - 0.6 / 250 s−1 One potential mechanism is cerebrospinal fluid cavitation occurring as a consequence of blast exposure or impact. Cavitation is the formation and expansion of vapor bubbles in a liquid when the hydrostatic pressure drops below the vapor pressure. The collapse of the bubble then causes a local shock wave that produces a compressive strain field onto nearby tissue, which may lead to damage. The bubble collapse can cause an even higher strain field than the strain field associated with the incident shock wave passage [123]. Previous studies have experimentally investigated cavitation generated by blast or blunt impact using SHPB [123], Shock Tube[153] [48], and FlyerPlate [24]. Multiple computational studies have also been performed to simulate the cavitation incident during TBI [93] [127]. The possible conditions leading to the low-pressure sites where cavitation occurs during TBI are a) in case of blast-induced TBI, the negative portion of the incident blast wave [93], b) skull-flexure [111], and c) negative wave reflection at the boundaries because of the acoustic impedance mismatch of tissue, skull, and CSF. It should also be noted that cavitation may not necessarily occur solely in the CSF. Cavitation has been observed within soft solid materials, which could indicate that it may happen within the brain tissue. 76 Despite an abundance of theoretical and experimental findings, there is as yet no unambiguous direct evidence of the occurrence of cavitation in a human brain during impact. Even if such a piece of evidence were to be found, it is hard to distinguish whether the injury results from cavitation or another effect pathologically. Still, the cavitation-induced damages is a reasonable hypothesis to explain the injuries seen in the contrecoup which is likely the site of cavitation. Multiple computational and experimental studies observed the pressure falls below the threshold in the contrecoup, [76, 79]. In the previous study [145], the author has computationally predicted the possible sites where cavitation occurs in the contrecoup during blast exposure. Accumulating experimental and computational evidence shows that cavitation could also occur within a sulcus during blast exposure [88] and impact [76]. Such an incident can cause a tauopathy seen in CTE. In that regard, the Mejia research group has been studying the intrasulcal cavitation and the response of sulci to cavitation-induced deformation using a simplified brain phantom. My contribu- tion to this effort focused on build a finite element model to simulate the deformation resulted from the cavitation occured in the intrasulcal region. The numerical model estimated that the maximum strain occurs at the depth of sulcus during the expansion of the cavitation bubble [76]. That might explain the occurrence of tauopathy at the depth of sulci. As additional background, an isolated single sulcus with varying sulcal depth and radius compressed by a rigid indenter was simulated to investigate high deformation regions within sulcal cavity in the work of Braun et. al [16]. Although the magnitude of maximum strain varies, depending on the geometrical parameters, the maximum strain was similarly observed at the depth of the sulcus in all cases. Regardless of the injury mechanism and source of impact, the brain soft tissue experiences a mechanical deformation, which initiates a cascade of neurotraumatic consequences following TBI. High-strain rate deformation is sufficient alone to lead to tau mislocalization and phosphorylation- dependent synaptic dysfunction that is the hallmark of CTE [16]. Current understanding of injury mechanisms of TBI across all spatial scales is still incomplete. Thus, there is a need for a multidisciplinary effort to understand how physical forces lead to injuries and affect the neurobiological features across the multiple scales of the brain. Following neuropathological 77 findings, it can be suggested that the pathological hallmark of CTE is localized mainly in the region from the depth of the sulcus spreading towards gray-white matter interfaces and the perivascular regions. The computational study presented in the rest of the chapter was developed to investigate the response of the brain tissue phantom to cavitation, observed in the experimental study [76]. Sections 4.2 and 4.3 presents the computational treatment and results previously published in [76]. 4.2 Computational model of cavitation-induced deformation The intrasulcal cavitation during the experiments grossly expand the sulci (cf. figure 4 on page 6 in [76]). The experimental results indicated that the onset for cavitation is consistent with the theoretical value of 𝑃𝑣 ≈ −98.19 kPa (gauge), calculated for water at 20◦C and at the local atmospheric pressure (101.32 MPa). The experimental onset for cavitation was determined based on the presence of the first nucleated bubbles in the field of view. Cavitation has been theorized to be an injury mechanism - often attributed to the high focal pressures after the bubbles collapse [17]. However, these experimental results also demonstrate the gyral deformation from bubble expansion within the sulci. From a mechanistic point of view, sulci could be considered as “notches” in a solid material. As such, their typical response to the effect of separating their walls from each other would be to concentrate strain / stress at their ends (sulcal depths). In order to semi-quantitatively evaluate the material response in the depths of the sulci, a series of numerical simulations of gyral response to bubble formation is performed. Experimental observations as presented in [76] led to the creation of an elastic two-dimensional (2D) computational model to simulate the response of the brain tissue phantom to the cavitation present in the blunt-impact experiments. To this aim, a numerical simulation of a single sulci was developed to examine the resultant tissue strain from intrasulcal bubble expansion. We assumed a plane strain and stationary state to simulate the deformation pattern resulted from the formation of cavitation bubbles. Figure 4.2 shows the domain of the computational model. Let X denote the initial position of a point in the given domain. The application of surface tractions will deform this material point X 78 to a new location x in a deformed configuration. The displacement of this point is defined by u ≡ x − X, (4.2.1) where x and X are vectors showing the initial and deformed location of the same material point respectively. Thus, the displacement is also a vector and explicitly given in u = 𝑢𝑥 (𝑥, 𝑦, 𝑧)ˆi + 𝑢𝑦 (𝑥, 𝑦, 𝑧)ˆj + 𝑢𝑧 (𝑥, 𝑦, 𝑧) ˆk (4.2.2) where 𝑢𝑥, 𝑢𝑦 and 𝑢𝑦 are the components of displacement vector u in the cartesian coordinate system with axes 𝑥, 𝑦, 𝑧 and unit vectors ˆi, ˆj, ˆk. A small strain treatment is appropriate when the displacements 𝑢𝑥, 𝑢𝑦, and 𝑢𝑧 are small. Thus, there is no distinction between the initial and deformed configuration, i.e., x ≈ X. In linear stress-strain analysis, the deformation is measured by a small strain tensor 𝝐 associated with the displacement field (4.2.2) that is given by 𝝐 = 1 2 (∇u + ∇uT) = , (4.2.3) 𝜖𝑥𝑥 𝜖𝑦𝑥 𝜖𝑧𝑥 𝜖𝑥𝑦 𝜖𝑦𝑦 𝜖𝑧𝑦 𝜖𝑥𝑧 𝜖𝑦𝑧 𝜖𝑧𝑧                     where ∇u denotes the displacement gradient. Considering a plane strain assumption where the displacement in the out-of-plane is equal to zero, the displacement vector (4.2.2) reduces to u = 𝑢𝑥 (𝑥, 𝑦)ˆi + 𝑢𝑦 (𝑥, 𝑦)ˆj + 0 ˆk. (4.2.4) Using (4.2.4) it follows that the strain components out of the plane are equal to zero i.e., 𝜖𝑧𝑧 = 𝜖𝑥𝑧 = 𝜖𝑦𝑧 = 0. Thus, the small strain tensor (4.2.3) 𝝐 = . (4.2.5)           𝜖𝑥𝑥 𝜖𝑦𝑥 𝜖𝑥𝑦 0 𝜖𝑦𝑦 0 0 0 0           At this point, a constitutive model to relate the deformation through the small-strain tensor (4.2.5) and the stress should be specified to conduct the analysis. The inner and outer layers are considered as isotropic linear elastic material. In the linear elasticity, the stress tensor 𝝈 is given by 𝝈 = C : 𝝐 ⇔ 𝝈ij = Cijkl : 𝝐kl, (4.2.6) 79 where C is the fourth-order elastic stiffness tensor and : represents the double contraction operation. For an isotropic linear elastic material, two material constants are enough to represent the elasticity tensor C. In this work, the bulk modulus 𝐾 and the shear modulus 𝜇 are employed. The fourth order stiffness tensor then has components Cijkl = (cid:18) 𝐾 − (cid:19) 𝜇 2 3 𝛿𝑖 𝑗 𝛿𝑘𝑙 + 𝜇(𝛿𝑖𝑘 𝛿 𝑗𝑙 + 𝛿𝑖𝑙 𝛿 𝑗 𝑘 ), (4.2.7) where 𝛿𝑖 𝑗 is the Kronecker delta (𝛿𝑖 𝑗 = 1 if 𝑖 = 𝑗 and 0 otherwise). One can obtain the other material parameters using the well-known conversion formulas. Due to the small-strain assumption, one does not make a distinction between all various stress tensors in finite deformation. Considering the plane strain assumption, the stress tensor (4.2.6) is given in the matrix form 𝝈 = 𝜎𝑥𝑥 𝜎𝑥𝑦 𝜎𝑦𝑥 𝜎𝑦𝑦 0 0 0 0 𝜎𝑧𝑧                     , (4.2.8) where the components of stress tensor is explicitly expressed as 𝜎𝑥𝑥 = (cid:0)𝐾 + 4𝜇 3 𝜎𝑥𝑦 = 2𝜇𝜖𝑥𝑦, 𝜎𝑦𝑦 = (cid:0)𝐾 − 𝜎𝑧𝑧 = (cid:0)𝐾 − 2𝜇 3 2𝜇 3 (cid:1)𝜖𝑥𝑥 + (cid:0)𝐾 − (cid:1)𝜖𝑥𝑥 + (cid:0)𝐾 + 2𝜇 3 4𝜇 3 (cid:1)𝜖𝑦𝑦, (cid:1)𝜖𝑦𝑦, (cid:1) (𝜖𝑥𝑥 + 𝜖𝑦𝑦). One can also notice that 𝜎𝑧𝑧 ≠ 0 despite the fact that 𝜖𝑧𝑧 = 0. The balance of linear momentum is given by div𝝈 + 𝜌b = 𝜌 𝜕u2 𝜕𝑡2 (4.2.9) (4.2.10) where 𝜌, and b denote the stress tensor, density, body forces respectively. Body forces b are regarded negligible, whereupon (4.4.1) for an equilibrium state reduced to div𝝈 = 0. 80 (4.2.11) The computational domain is 40 mm long by 29.32 mm wide and discretized by nearly 14,300 quadrilateral elements. The initial sulcus gap was 1.1 mm, as shown in the blue shaded area in figure 4.2. The sulcus gap is also the location where the cavitation bubbles were observed from the high- speed imaging. The boundary consists of five parts. The most important is the intrasulcal boundary 𝜕Ω𝐼 which separates the gray matter in the sulcus with the sulcal cavity containing the CSF. In other words, the CSF itself is not part of the computational domain. The other four parts of the boundary are top, bottom, left,and right portions labelled 𝜕Ω𝑇 , 𝜕Ω𝐵, 𝜕Ω𝐿, 𝜕Ω𝑅 respectively, as shown in Figure 4.4. The presented computational model is developed using COMSOL Multiphysics ® v.5.6 Solid Mechanics Module. The solid portion which is actually modelled consists of two subdomains: an outer layer of gray matter and an inner layer of white matter, both of which were modeled using an isotropic linear elastic material model. The material properties were matched to that of the phantom, specifically the density and high-frequency (25 Hz) shear modulus of the PAA gelatin [142]. The other elastic properties (Bulk and Young’s modulus) were computed by well-known elastic conversion formulas given K = 2𝜇(1 + 𝜗) 3(1 − 2𝜗) , E = 2𝜇(1 + 𝜗) (4.2.12) where 𝜇, 𝜗, and E denote shear modulus, Poisson’s ratio, and Young’s modulus of the material. With the assumption of Poisson’s ratio of 𝜗 = 0.49, the material properties are tabulated in table 4.2. Table 4.2 The material parameters for Layer 1 and Layer 2 used in the phantom. Layer 1 (White) 2 (Gray) Density (kg/m3) 1056 1037 Shear Modulus (kPa) 4.19 7.27 Poisson’s Ratio 0.49 0.49 Bulk Modulus (kPa) 208.45 361.47 Young’s Modulus (kPa) 12.48 21.66 The effect of the initiation of cavitation-induced bubbles on the surrounding phantom was characterized by a pressure boundary condition on 𝜕Ω𝐼. The high-speed video images demonstrate that the range of total displacement occurring in the phantom adjacent to the cavitation bubble is 81 between 2.30 mm and 2.69 mm. A preliminary study was first performed to compute the pressure magnitude that can lead to the maximum displacement between two sulcal walls measured from the experiment. In the computational model, the boundary load F𝐴 was then modeled by applying the scalar pressure magnitude 𝑝 to the normal to the intrasulcal boundary i.e., F𝐴 = −𝑝n on 𝜕Ω𝐼 illustrated in figure 4.4. The magnitude of pressure 𝑝 applied to the sulcal cavity was 11.06 kPa. This resulted in 2.36 mm maximum displacement between the sulcal walls which is in the range of displacement observed in the experiment. In 3-D elasticity, three boundary conditions are required at each boundary point. Considering the plane strain model problem at hand, the number of boundary conditions reduces to two conditions. These two boundary conditions were defined to all outer boundaries of the domain except for the intrasulcal area represented by blue in figure 4.2. These are 1) the displacement is zero in the normal direction to the boundary u · ˆn = 0, (4.2.13) where ˆn is the outward unit normal, and 2) no shear traction, which is formally expressed as 𝝈 ˆn − ( ˆn · 𝝈 ˆn) ˆn = 0 on 𝜕ΩT, 𝜕ΩB.𝜕ΩR.𝜕ΩL. (4.2.14) Together the two boundary conditions (4.2.13) and (4.2.14) are sometimes depicted as a "roller type" boundary condition. Boundary condition (4.2.14) is a natural boundary condition, but (4.2.13) is a kinematic boundary condition. For this reason, (4.2.13) is explicitly enforced in numerical treatments, whereas (4.2.14) may not need to be stated depending on the computational routine. At the interface between white and gray matter, a displacement continuity condition was specified over the interface where layer 1 and layer 2 match. The presented model also includes traction continuity, which is a natural boundary condition because this condition naturally emanates from the weak form of (4.2.11). 4.3 Examination of the equilibrium response to pressure loading The numerical model allowed us to detect the locations where heightened mechanical strain is induced by sulcal expansion arising from cavitation-induced blunt trauma. The contours of strain 82 Figure 4.2 Geometrical illustration of 2D plane-strain single sulcus simulation domain [76]. are shown in figure 4.3 to illustrate the response of the brain’s soft tissue to cavitation-induced sulcal expansion. The results estimated a high tensile strain concentration at the sulcal depth due to the expansion. The maximum strain 𝜖𝑥𝑥 is measured 0.41% at the tissue / CSF interface specifically right at the deepest point of the sulcus suggesting that this is a point of stress concentration. The strain then reduces toward the bulk of the tissue. Based on figure 4.5, the strain decays along the 𝑦−axis when traversing from the tip of the sulcus toward the bulk of the tissue. Then, the strain 𝜖𝑥𝑥 transitions from tensile to compressive at 𝑦 ≈ 1.68 mm, and keeps decaying at a lower rate until it reaches the layer 1 / layer 2 interface at 𝑦 = 4.15 mm. At that point, the strain 𝜖𝑥𝑥 exhibits a slope discontinuity and reversal, relaxing asymptotically toward zero strain. Hence, this interface is also the location at which the strain attains its minimum value (−0.17), which is the maximum magnitude of compressive strain. As such, it can be considered a stress-concentrating feature, though milder than the stress-concentration at the tip of the sulcus. A qualitative comparison between our predicted contours of tensile strains secondary to bubble expansion and the tau pathology found at the depths of the sulci have multiple similarities as has 83 CavitationRegionYXLayer 2PAA 12 wt.%Layer 1PAA 10 wt.%LTBIR Figure 4.3 Computed normal strain 𝜖𝑥𝑥 around the sulcus due to deformation induced by expanding cavitation bubbles. The computational analysis solves an equilibrium boundary value problem that is driven by the pressure 𝑝 =11.06 kPa within the sulcal cavity. Red represents tensile strain, while blue represents compressive strain [76]. been reported in neuropathologic investigations of human subjects [40, 84]. Firstly, high tensile strains are located in the same location and have the same length scale ∼ 1𝑚𝑚 that is roughly symmetric between the two sides of the sulci and is of similar positioning as tauopathy observed in early CTE. The highest strain resulting from the expansion is measured as 0.41 from the location 84 Strain0.410.300.200.100-0.10-0.17xyxy0.00.554.151.1 mm3.6 mm of the deepest point of the sulcus. The compressive maximum strain in x-direction 𝜖𝑥𝑥 is 0.17 and measured around the endpoint of the semi-circular of the sulcus on either side, as shown with blue in Figure 4.3. The compressive strain seen in the simulation results are symmetric in nature and extend distally outwards from the deepest point of the sulcus on either side. The maximum compressive strain is observed at a ±75.6◦ angle off of the innermost portion of the sulcus. Figure 4.4 indicates the distribution of strain along the semi-circular CSF-brain interface at the sulcal depth. To characterize the distribution of strain along the semi-circular interface , a polar coordinate system with radial origin at the center of the semi-circular CSF-brain interface is utilized and 𝜃 = 0◦ at the middle point of the arc (deepest point of sulcus). Note that, for convenience, the angle as positive in the clockwise direction is defined. Given that the sulcus is 1.1 mm wide, the radius of the CSF-brain interface (sulcal depth) is 𝑟 = 0.55 mm. Based on this coordinate system, we recovered the strain as a function of angle along the CSF-brain interface in figure 4.4. From the figure, the maximum strain (0.41) is seen at 𝜃 = 0◦. The strain then reduces along the CSF-brain interface on both sides of the maximum, reaching a minimum strain value of −0.17 (maximum compressive strain) at ±75.6◦. Figure 4.4 Computed distribution of the normal strain component in 𝜖𝑥𝑥 along the semi-circular CSF-brain interface at the sulcal depth (r = 0.55 mm, −90◦ ≤ 𝜃 ≤ 90◦) [76]. 85 -0.2-0.100.10.20.30.40.5StrainCSF-brain interfaceSulcus’deepest pointMaximum tensile strainMaximum compressive strainSulcus’deepest pointLayer 1Layer 2 Figure 4.5 Computed distribution of the normal strain component in 𝜖𝑥𝑥 along the y-axis (probing line) [76]. Previous computational models have predominantly centered their attention on the localization of strain and strain rate at the depths of the sulci as a mechanism leading to a physical injury, [16]. While the physical injury mechanisms are mostly linked to strain and strain rates, comprehending the accumulation of stress is also important to understand these physical injury mechanisms. To this end, von Mises stress 𝜎mises is utilized i.e., 𝜎mises = 2s : s where s is the deviatoric 𝛿𝑖 𝑗 𝜎𝑘 𝑘 . The contour of figure 4.6 shows von Mises stress resulting from the stress 𝑠𝑖 𝑗 = 𝜎𝑖 𝑗 − 1 3 √︃ 3 deformation triggered by the expansion of a cavitation bubble within the intrasulcal gap. Likewise, the results clearly indicate that the stress field is localized particularly within the sulcal arc and its adjacent walls. The maximum von Mises stress is recorded as 5.3 kPa along the semi-circular interface between the cerebrospinal fluid (CSF) and gray matter, making this region a site of stress concentration. Subsequently, the von Mises stress levels diminish as one moves deeper into the bulk of the brain tissue. 86 246810y-axis (mm) (probing line)-0.100.10.20.30.40.5StrainLayer 2Layer 2xy0.554.155.0010.0Layer 1The Layer 1-2 Interface00The CSF-brain Interface4.150.55Slope discontinuityLayer 1 Figure 4.6 Von-Mises stress pattern within the single sulcus domain. 87 1.8245.3Von Mises (kPA) 4.4 Examination of the response of sulci to transient loading In the preceding section, a static boundary condition is implemented to examine the response of a single sulcus to the deformation induced by the expansion of a cavitation bubble within the intrasulcal region. However, it is essential to acknowledge that the formation and expansion of cavitation bubbles represent a time-dependent phenomenon, as demonstrated in Figure 4 of the work by [76]. To address the temporal response of the material, a 2D plane strain computational model is developed in the context of linear elastodynamics. For this purpose, a time-dependent loading profile is adopted to simulate the expansion of the intrasulcal bubble. In the absence of the body forces, the balance of linear momentum (4.2.10) is reduced to div𝝈 = 𝜌 𝜕u2 𝜕𝑡2 . (4.4.1) Following the framework presented in section 4.2, the geometry of reference configuration, mesh structure, elastic properties and the boundary conditions outside the sulcus remain the same in this fully-dynamic computational model. The primary modification lies in the specification of the time-dependent pressure boundary condition on the intrasulcal cavity boundary 𝜕Ω𝐼. Zero displacement initial conditions are taken throughout the domain. The pressure loading is in the form of a ramp and hold function as shown in the figure 4.7. The slope of this linearly increasing portion of the pressure profile on this boundary is obtained from experimental data spanning from 0.765 ms to 1.761 ms, as presented in Figure 5 in [76]. The rate of pressure increase is around 5.5𝑥103 kPa/ms. It is significant to note that the experimental data was obtained from a pressure sensor located in the counter-coup region, and it is assumed that the characteristics of the pressure profile developing within the intrasulcal region are akin to those in the counter-coup region. A linear elastic constitutive model is again utilized to model the gray and white matter, and as such this analysis does not treat time dependent effects within the constitutive response itself, such as viscoelastic or poroelastic phenomena. The simulation was conducted over the initial 7.5 ms subsequent to the initiation of cavitation bubble expansion. 88 Figure 4.7 Pressure loading profile applied to the external boundary 𝜕Ω𝐼. In what follows we specifically examine the quantity known as the equivalent deviatoric strain, denoted as 𝜖𝑑𝑒𝑣 which is 𝜖𝑑𝑒𝑣 = 2 3 √︃ 𝑥𝑥 − 𝜖𝑥𝑥𝜖𝑦𝑦 + 𝜖 2 𝜖 2 𝑦𝑦 + 3𝜖 2 𝑥𝑦. The equivalent deviatoric strain combines multiple strain tensor components into a scalar value. Hence, it facilitates examine the regions experiencing a strain in various directions. The contours of figure 4.8 illustrate the temporal and spatial distribution of 𝜖𝑑𝑒𝑣 the deformation within the sulcus which specifically accumulated in the regions of gray matter and the gray-white matter interface. The highest values of 𝜖𝑑𝑒𝑣 are recorded around the semi-circular interface between cerebrospinal fluid (CSF) and gray matter at the depth of the sulcus throughout the deformation history. Furthermore, the sulcal walls adjacent to the sulcal arc exhibit relatively elevated deviatoric strain during the deformation. This is represented in Figure 4.8, which also illustrates the deforma- tions undergone by the single sulcus in response to the pressure ramp. The black boundary lines within Figure 4.8 delineate the reference configuration prior to any deformation. In Figure 4.9, the change in 𝜖𝑑𝑒𝑣 over time in the deepest point of the sulcus, which corresponds to 𝜃 = 0 shown in the inset of figure 4.4. The 𝜖𝑑𝑒𝑣 experiences an increase within the first 2.42 89 .Level of static loading Figure 4.8 Time-lapse images of 𝜖𝑑𝑒𝑣 at nine different time step. ms of the simulations, reaching a local maximum value of 0.47. During this time interval, both the sulcal walls and the sulcal arc expand, as evident in the top row of figure 4.8. It is followed by the propagation of the deformation along the sulcal walls moving away from the sulcal arc, 90 t = 1.0 mst = 1.5 mst = 2.0 mst = 3.0 mst = 3.5 mst = 4.0 mst = 4.5 mst = 6.0 mst = 7.5 ms0.350.30.250.20.150.10.050Equivalent Deviatoric Strain which gradually returns to its equilibrium position simultaneously. At time instance 𝑡 = 3.82 ms and 𝑡 = 4.69 ms, the sulcal arc region exhibits the lowest levels of deviatoric strain. However, the sulcal arc region starts expanding again which can be explicitly seen in figure 4.8 at 𝑡 = 6.0 ms. This resurgence in expansion corresponds to another surge in equivalent deviatoric strain, reaching a maximum value of 0.52 at 𝑡 = 6.67 ms. Figure 4.9 Computed time history of equivalent deviatoric strain observed at the deepest point of the sulcus. One potential underlying physical mechanism associated with CTE involves the localization of strain and strain rate within the depths of the sulci during brain motion following an impact. Hence, our focus now shifts toward an analysis of the temporal and spatial distribution of the strain rate induced by the deformation resulting from the expansion of the cavitation bubble. The rate of strain tensor is considered to address the possible effect of the strain rate. The rate of strain tensor is the symmetric part of the spatial gradient of velocity field i.e., 𝑑𝑖 𝑗 = 1 2 (𝑙𝑖 𝑗 + 𝑙 𝑗𝑖) where the spatial gradient of velocity field 𝑙𝑖 𝑗 is given by 𝑙𝑖 𝑗 = 𝑣𝑖, 𝑗 . Figure 4.10 demonstrates the variations in the normal components of 𝑑𝑖 𝑗 , as observed at the deepest point of the sulcus. Both 𝑑𝑥𝑥 and 𝑑𝑦𝑦 exhibit 91 a fluctuating pattern over the course of the simulation and attains local maxima at time instances of 1.10, 3.66, and 5.69 ms. The maximum strain rate is measured at 429 (1/s) in the x-direction and -456 (1/s) in the y-direction at 3.65 ms. Figure 4.10 Computed time history of the normal components 𝑑𝑥𝑥 (solid - black) and 𝑑𝑦𝑦 (dashed- red) at the deepest point of the sulcus. For the aforementioned time instances, the spatial distribution of strain rate is illustrated in the contours depicted in Figure 4.11, which provides insight into the spatial representation of the strain rate at the time when the maximum strain rate is observed. It should also be noted that figure 4.11 illustrates the deformations as well and the black boundary lines delineate the reference configuration prior to any deformation for the given time instance. The gray matter experiences relatively higher strain rate comparing to the white matter during the computed time histories, specifically localized around the sulcal arc and the sulcal walls. The results reveal that the deformation generated by the expansion of the cavitation bubble at the intrasulcal region leads to highly localized strain and strain rate fields at the deepest of the sulcus. From a mechanistic 92 (a)(b)(c)(e)(f)(d)dd1.103.655.66 Figure 4.11 Time-lapse images of normal components 𝑑𝑥𝑥 (top row) and 𝑑𝑦𝑦 (bottom row) directions at 1.10, 3.65 and 5.66 ms (from left to right). point of view, due to the anatomical structure of the sulcus and the loading conditions, the u-shaped semicircular notched edge analogy might be adopted to explain the highly localized strain pattern in the depth of the sulcus. Previous studies [99] show that there is a relation between the geometrical properties of the notch and the strain / stress concentration factor. The ratio of the depth of the sulcus to the radius of the semi-circular end is linearly proportional to the concentration factor. However, it has not been investigated if this relationship holds for bubble-induced strain in sulci – as found in this study. For an isotropic linear elastic material there are two characteristic wave speeds in the bulk. 93 0x102321-4-3-2-1StrainRate (1/s)(a)(b)(c)(d)(f)(e) The first is the shear wave speed 𝑐𝑠 = √︃ 𝜇 𝜌 which is characterized by the perpendicular particle displacement to the wave propagation direction. The second is the dilatational wave speed 𝑐𝑑 = √︃ 3𝐾+4𝜇 3𝜌 which is characterized by the particle displacement in the direction of wave’s propagation [87]. For the properties given in Table 4.2 these wave speeds are shown in Table 4.3. With the aid of these elastic wave properties, the distance between the gray-white matter interface to side boundaries (e.g., 𝜕Ω𝐿, 𝜕Ω𝑅) and top boundary 𝜕Ω𝑇 , and the thickness of gray matter which are 14.11 mm, 17.91 mm and 3.6 mm respectively, one can compute how far waves travel during the total simulation time. This in turn enables the determination of specific times at which first interactions take place. Namely, dilatational waves originating at the external boundary 𝜕Ω𝐼 then travel through the gray matter, reaching the white matter interface after 0.19 ms. The required time for the dilatational waves to reach the lateral boundaries and reflect back to the sulcal cavity boundary 𝜕Ω𝐼 is around 2.36 ms. Hence, the dilatational waves can travel multiple times back and forth between the lateral boundaries (e.g., 𝜕Ω𝑅, 𝜕Ω𝐿), and external boundary 𝜕Ω𝐼. The dilatation waves reflected from the lateral boundaries reach the gray and white matter interface at 2.17 ms, 4.53 and 6.89 ms and the sulcal cavity at 2.36 ms, 4.72 ms and 7.08 ms during the simulation time. Similarly, there is a sufficient time for the dilatational waves to reach the top boundary 𝜕Ω𝑇 and reflect back to sulcal cavity boundary 𝜕Ω𝐼. Due to the longer distance between these two boundaries, the required time for the dilatational wave to travel back and forth is slightly different, which is around 2.89 ms. Therefore, the reflected dilatational waves reach to the sulcal cavity at 2.89 ms and 5.79 ms during the simulation time. Similarly, shear waves originating from the sulcal cavity 𝜕Ω𝐼 reaches the white matter interface after traveling through the gray matter. Due to the relatively lower wave propagation velocity, shear waves reach the interface around 1.36 ms, then continue to propagate through the far boundaries. However, the shear waves interacting with the interface can propagate around 12.21 mm within the simulation time considered. In other words, the shear waves do not have sufficient time to regain the gray-white matter interface and sulcal cavity boundary 𝜕Ω𝐼. Hence, shear wave interaction 94 around the interface and the sulcal depth was not observed in this model. The impedance of both subdomain associated with the dilatational 𝑍𝑑 and shear 𝑍𝑠 wave are computed to understand how the waves behave specifically at the gray-white matter interface, i.e., 𝑍 = 𝜌𝑐 where 𝑐{·} denotes the shear 𝑐𝑠 and dilatational 𝑐𝑑 wave velocity. The specific impedance of gray and white matter are also given in Table 4.3. The impedance of gray matter associated with the dilatational and shear wave velocity are approximately 26% greater than that of white matter. This impedance difference results in reflecting of certain part of the incident wave while the rest part is transmitted through the interface. The reflection coefficient 𝑅 at the gray-white matter interface was computed -0.13, i.e., 𝑅 = 𝑍𝑤−𝑍𝑔 𝑍𝑔+𝑍𝑤 where 𝑍𝑔 and 𝑍𝑤 are the specific impedance of gray and white matter respectively. The transmission coefficient 𝑇 at the gray-white matter interface was computed around 0.86, i.e., 𝑇 = 2𝑍𝑤 𝑍𝑔+𝑍𝑤 . These results indicates 13% of an incident wave that is normal to the interface is reflected while 86% of the wave is transmitted. The sign of the reflection and transmission coefficients shows whether the phase of reflection and transmission waves are in the same phase with the incident wave. In other words, the reflected wave is 180◦ out phase with the incident wave whereas the transmitted wave is in the same phase with the incident wave at the interface. Table 4.3 The elastic wave properties for Layer 1 and Layer 2 used in the phantom. Layer 1 (White) 2 (Gray) 𝑐𝑑 (m/s) 14.23 18.91 𝑐𝑠 (m/s) 1.99 2.64 𝑍𝑑 (kg/m2s) 15026.88 19609.67 𝑍𝑠 (kg/m2s) 2103.48 2745.72 The above commentary on the reflection and transmission of the waves is based on impedance mismatch between the layers and the incident waves that are normal to the interface. In such solid- solid interfaces, mode conversion also occurs when an incident wave encounters an interface with an oblique angle. Depending on the angle of the incident wave, the impedance of the layers, and the type of the wave, mode conversion takes place from dilatational wave to shear wave or vice-versa. The angle of refracted dilatational and shear waves can be computed utilizing Snell’s law. Such interactions can also give rise to localized interface waves as the angle of the incident wave becomes 95 relatively shallower. In other words, this leads to a dilatational wave that travels horizontally along the interface while all incident waves are reflected. The critical angle of incident wave, which causes an interface wave, can be similarly computed using Snell’s law for the case where the angle of refraction is 90◦, i.e., sin(𝜃1) 𝑐𝑑 (1) where 𝑐𝑑 (1) and 𝑐𝑑 (2) represent the dilatational wave velocity = sin(𝜃2) 𝑐𝑑 (2) of layer 1 and 2 while 𝜃1 denote the angle of incident wave, and the angle of refracted wave 𝜃2 = 90◦. For the considered bilayer simplified sulcus model, the critical angle 𝜃𝑐𝑟 is 49◦ with respect to the normal to the interface. In other words, any incident wave traveling from the white matter to gray matter with an angle 𝜃1 ≤ 49◦ could potentially lead to an interface wave. It is important to note that the computational model assumes the reference configuration of the sulcus to be stress-free. However, the biological growth and remodeling alters the mechanical state of the soft tissues and leads to residual stress and strain fields. These residual stress fields are inhomogeneous, and anisotropic (including both compression and tension components due to the zero-traction condition), thereby leading to more complex mechanical behavior compared to stress-free conditions. Furthermore, these residual stress fields have the potential to significantly influence the baseline response of soft tissue to external forces, as presented in Chapter 2. During the process of cortical growth and gyrification, the brain’s soft tissue experiences large strains as a consequence of inhomogeneous growth. In the context of brain injury dynamics, it is imperative to consider the effects of residual stress fields for a more accurate understanding of how the brain’s soft tissue responds to various external loadings. Hence, this study proceeds with the examination of growth-induced residual stress fields that emerge within sulcus-gyrus formation, aiming to gain more insights into the reference mechanical state of the brain’s soft tissue. 96 CHAPTER 5 THE GROWTH-INDUCED RESIDUAL STRESS FIELDS RESULTING FROM SULCUS-GYRUS FORMATION IN THE BRAIN 5.1 Preliminaries and kinematics of finite volumetric growth The kinematics of finite volumetric growth takes into account the combined effect of growth and elastic deformation, which was originally proposed in [109]. Let us assume a stress-free unloaded reference configuration denoted by B0 in which X is a position vector prior to any deformation. Growth leads to the change in shape and size motivated by various biological mechanisms, that is regarded as given for the present treatment. The growth is kinematically described by a tensor valued function Fg mapping the reference configuration B0 into a virtual intermediate configuration B𝜉. The growth field 𝑣 does not need to be continuous across the body that leads to discontinuity within the body. Hence, the intermediate configuration necessarily undergoes an elastic deformation to restore compatibility subsequently. The elastic deformation tensor maps the virtual intermediate configuration B𝜉 to a deformed configuration B. Equivalently, the growth followed by elastic deformation is mapped by a total deformation gradient F which is given in (1.3.1). Following [109], the total deformation gradient F is decomposed into the growth tensor and the elastic deformation tensor, as given in (1.3.2) The energetic framework of hyperelastic materials is represented by a stored energy density 𝑊 as a function of Fe. The strain energy density depends on Fe only via the right Cauchy-Green tensor Ce due to the frame invariance of strain energy density i.e., 𝑊 = 𝑊 (Ce, X) where Ce is the right Cauchy-Green tensor associated with the elastic accommodation Fe as expressed in (1.3.5). In the context of morphogenesis, the growth process occurs on a relatively long time scale ranging from days to months which is larger than any other time scales (e.g., relaxation time scale of soft tissue). Hence the process can be considered quasi-static. With the absence of the body forces, the equilibrium equation thereby reduces to div T = 0, 97 (5.1.1) where div is the divergence operation with respect to the deformed location x and T is the Cauchy stress. The elastic deformation accommodating the growth introduces internal mechanical stress within the body, which preserve equilibrium by satisfying (5.1.1) without any surface tractions. Confining our attention to compressible isotropic hyperelastic materials, the strain energy function 𝑊 depends upon Ce only through its principal scalar invariants 𝐼 𝑒 and the third invariant 𝐼 𝑒 e)] 3 = detCe = 𝐽𝑒 denotes the volumetric change due to the elastic deformation. The strain energy function for a nearly incompressible isotropic hyperelastic material is expressed in the form of 𝑊 = 𝑊DEV(𝐼 𝑒 1 ) + 𝑊VOL(𝐽𝑒) where 𝑊VOL represents the pressure-volume 2 [(tr Ce)2 − trC2 1 = tr Ce, 𝐼 𝑒 2 = 1 , 𝐼 𝑒 3 , 𝐼 𝑒 2 response. The Cauchy stress tensor is T = 2 𝐽𝑒 Fe 𝜕𝑊 𝜕Ce FT e = 2 𝐽𝑒 (cid:18) 𝜕𝑊 𝜕𝐼 𝑒 1 + (cid:19) 𝐼 𝑒 1 𝜕𝑊 𝜕𝐼 𝑒 2 Be + (cid:19) 2 𝐽𝑒 (cid:18) 𝜕𝑊 𝜕𝐼 𝑒 2 B2 e + 𝜕𝑊 𝜕𝐽𝑒 I (5.1.2) where Be = FeFT e is the left Cauchy-Green deformation tensor. Equation (5.1.2) describes growth- induced mechanical stresses building up in an isotropic compressible hyperelastic material [109, 50]. Figure 5.1 Schematic representation of kinematics of finite volumetric growth adapted to the cortical growth and folding process (top). Idealized outer layer-inner core system of the computational domain with assigned geometrical parameters, and details of the boundary conditions (bottom). The hypothesis under consideration is that brain morphogenesis is driven by mechanical insta- bilities in a given bilayer configuration that arise from the growing mismatch between an expanding 98 freesymmetric symmetric clampedΩXcΩXL∂ΩXT∂ΩXI∂ΩXR∂ΩXB∂ΩXsΩXCurrent State(Residually-stressed)Reference(stress-free)IntermediateStateOuter LayerInner CoreYXHLFFeFgHscTangentialExpansion outer layer and the relatively less-growing inner core. As depicted in figure 5.1, this simply pro- poses that the outer layer grows more than the inner core in the tangential direction. The differential growth of the layers leads to incompatible deformation fields as depicted in the intermediate con- figuration B𝜉. Elastic deformation Fe that restores the compatibility of the body subsequent to the growth leads to a mismatch strain 𝜀𝑀 in the outer layer 𝜀𝑀 = (𝐿𝑜 − 𝐿𝑖)/𝐿𝑜, (5.1.3) where 𝐿𝑜 and 𝐿𝑖 is the arc length of outer layer and inner core in the intermediate configuration B𝜉. Depending on the material and geometrical properties, the mismatched strain 𝜀𝑀 could cause various patterns (e.g., folding, wrinkling, crease etc.) in the surface of the current configuration B. Considering our particular interest in understanding the internal mechanical state of a gyrus- sulcus formations, the inner and outer layer represent the subcortex (white) and cortex (gray matter), respectively. The volume of the cortex and subcortex increase over time for various reasons (e.g., cell growth, proliferation, accumulation, etc.). However, the underlying biological mechanism is beyond the scope of this work. We proceed by building a computational model to find out the internal mechanical states of the soft brain tissue in response to growth and gyrification qualitatively. To this end, a portion of the brain is idealized as a rectangular block to capture the mechanics of gyrus-sulcus formation. Considering the block in its reference configuration occupying the region ΩX before any growth that is given by Cartesian coordinates 𝑋, 𝑌 , 𝑍 where 0 ≤ 𝑋 ≤ 𝐿, 0 ≤ 𝑌 ≤ 𝐻, −∞ ≤ 𝑍 ≤ ∞. (5.1.4) Here, 𝐿 and 𝐻 denote the length and total height of the domain, i.e., 𝐻 ≡ 𝐻𝑐 + 𝐻𝑠 where 𝐻𝑐 and 𝐻𝑠 represent the inital height of the cortex and subcortex as depicted in figure 5.1. Then, the reference configuration ΩX is partitioned into two sub-regions i.e., ΩX = ΩX𝑐 ∪ ΩX𝑠 including an subcortex ΩX𝑠 with 0 ≤ 𝑌 ≤ H𝑠 and cortex ΩX𝑐 with H𝑠 ≤ 𝑌 ≤ H𝑐. 99 5.2 The computational model of Growth-Induced Surface Instabilities In the presented treatment, each layer requires growth tensor Fg with the distinct growth factors to define differential growth. The growth tensor is expressed for both layers in the form of Fg = 𝛾𝑋E𝑋 ⊗ E0,𝑋 + 𝛾𝑌 E𝑌 ⊗ E0,𝑌 + 𝛾𝑍 E𝑍 ⊗ E0,𝑍 , (5.2.1) where {E0,𝑋, E0,𝑌 , E0,𝑍 } and {E𝑋, E𝑌 , E𝑍 } are the unit basis vectors in the reference B0 and inter- mediate B𝜉 configurations, while 𝛾𝑋, 𝛾𝑌 and 𝛾𝑍 represent the growth factors in the corresponding directions. Following [50], two basis can be identified at each point so that E𝑋 = E0,𝑋, E0,𝑌 = E𝑌 , E0,𝑍 = E𝑍 . The growth factor should be greater than the unity to define positive growth (expansion) in the associated direction. If the growth factor equals unity, any growth does not occur within the body. In case of the growth factors are less than unity, it leads to shrinkage (or atrophy). Following the framework described in section 2, an elastic deformation requires to ensure the compatibility of the layers subsequent to the growth. That enables a compatible displacement field in the body u. The total deformation gradient (1.3.1) is thus expressed F = I + 𝜕u 𝜕X . (5.2.2) 5.2.1 Variational treatment It is assumed that the resultant folded pattern always seeks the lowest potential energy state among all other kinematically admissible states. In the absence of body forces, equilibrium deformations are determined on the basis of choosing u that minimizes the energy formalism Π = Πstored − Πload. Here, Πload is the work due to the external loads and the total stored elastic energy Πstored is the integral of strain energy density per unit length 𝑊 (Fe) over the volume of reference configuration Πstored = ∫ ΩX 𝑊 (Fe) 𝐽𝑔 d ΩX, (5.2.3) where 𝐽𝑔 = det F𝑔 denotes the volume change following growth. For the bilayer system de- picted in figure 5.1, the total stored elastic energy (5.2.3) can be explicitly expressed by Πstored = ∫ Ωx𝑠 𝑊 (Fe) 𝐽𝑔,𝑐 d ΩX𝑐 where 𝐽𝑔,𝑠 and 𝐽𝑔,𝑐 denotes the volume change due 𝑊 (Fe) 𝐽𝑔,𝑠 d ΩX𝑠 + ∫ Ωx𝑐 100 to the growth in the associated subdomain. In the computational model, we seek the displacement field u minimizing the total stored energy Πstored which corresponds to the equilibrium state of each grown configuration u = arg min u∈𝑉 Πstored, (5.2.4) where 𝑉 is the function space consisting of first-order continuous Lagrange finite element functions that satisfy the boundary conditions. It should also be noted that the variation of (5.2.3) with respect to the displacement field yields to the weak form of linear momentum balance with the absence of the surface tractions. 5.2.2 Growth stipulation The presented work proceeds with a special case to determine the qualitative behavior of growth-induced residual fields. To this aim, an initially flat bilayer model as depicted in Figure 5.1 is employed. The bilayer model given in (5.1.4) is assumed to be infinitely thick and, thus, it imposes its in-plane strains. The initial thickness of the cortex 𝐻𝑐 is selected 2 mm, and the subcortex thickness is ten times as thick as the cortex. The length of the domain includes the half-wavelength of the fold with symmetric boundary conditions. It is therefore determined based on (1.2.1) for a given shear modulus contrast, and the initial cortex thickness 𝐻𝑐. Following the differential tangential growth hypothesis to describe the growth kinematics of the brain soft tissue, only tangential (in-plane) growth in the cortex is considered whereas the subcortex does not undergoes any growth in any direction. With the reduction from 3D to 2D, the same notation for Fg given in (5.2.1) is maintained to use after dropping z-component. The growth factors are specified as 𝛾𝑋 = (1 + 𝑔) and 𝛾𝑌 = 1 for the cortex where 𝑔 is a pure scalar that controls the amount of the growth in the tangential direction. Note that the parameter 𝑔 is strictly positive 𝑔 > 0. The growth tensor (5.2.1) thereby becomes Fg = (1 + 𝑔)E𝑋 ⊗ E𝑋 + E𝑌 ⊗ E𝑌 . The growth factors are 𝛾𝑋 = 𝛾𝑌 = 1 for the subcortex, which means that any growth does not occur. Due to the considered growth conditions, the cortex experience distinct change in size, resulting in a compressive mismatched strain in the cortex to maintain the compatibility of layers. Recalling the 101 definition given in (5.1.3), the mismatch strain 𝜀𝑀 turns into for the considered growth condition 𝜀𝑀 = 𝑔/(1 + 𝑔). (5.2.5) As the mismatch strain 𝜀𝑀 increases, the compressive stress continues to build up in the cortex. Once the mismatch strain exceeds a certain threshold, the cortex buckles out of the growth plane to release accumulating compressive stress, resulting in various surface instability. Similarly, the compatible displacement field is reduced to form of u = 𝑢𝑥 (𝑥, 𝑦)E𝑥 + 𝑢𝑦 (𝑥, 𝑦)E𝑦 for the considered model where out-of-plane displacement is zero. Here, the displacement field u, which minimizing the total storage energy (5.2.4), is limited with the length in the 𝑥-direction whereas the upper bound for the displacement field in the y-direction is determined by the parameter 𝐷, i.e, 0 ≤ 𝑢𝑥 ≤ 𝐿 and 0 ≤ 𝑢𝑦 ≤ 𝐻 + 𝐷. In the context of brain morphogenesis, the parameters 𝐷 can be considered the distance between the top surface of cortex and the inner surface of the skull. Although the skull is not required to trigger the gyrification [8], previous computational works suggest that the existence of the skull may affect the final morphology of the cortex by flattening the gyral crowns [125, 90]. While this approach facilitates the incorporation of potential influences stemming from the skull, the skull constraint is neglected until for a later consideration. 5.2.3 Boundary conditions External loads (e.g., the pressure generated by the surrounding intracranial structures such as cerebrospinal fluid (CSF)) will initially be disregarded in the computational model because such an external loading act as a mechanical perturbation that can potentially affect the region where the folding initiates and lead to a different folded morphology. Along with external loads, any anatomical and mesh imperfection, which could act as initial perturbation, were not introduced in the model either. The top surface of outer layer 𝜕ΩT X is thus traction free. The boundary conditions on the right 𝜕ΩR X and left walls 𝜕ΩL X of the domain are: i) the displacement is zero in the normal direction to the boundary to preserve symmetry u · ˆn = 0 on 𝜕ΩR X , 𝜕ΩL X (5.2.6) 102 where ˆn is the outward unit normal, and ii) no shear traction, which is formally expressed as T ˆn − ( ˆn · T ˆn) ˆn = 0. (5.2.7) The bottom surface is held fixed (i.e., u = 0 on 𝜕ΩB X) to provide a good approximation of a half- space for the inner core. The nature of implementation does not explicitly require enforcement of zero traction boundary condition since this is a natural condition associated with the variational procedure. The boundary conditions are also illustrated in figure 5.1. 5.2.4 Interface conditions We restrict consideration to the case where the cortex and subcortex are perfectly bonded. The displacement field is to be continuous at the interface u(𝐻+ 𝑠 ) = u(𝐻− 𝑠 ) where 𝐻+ 𝑠 and 𝐻− 𝑠 represent the outer layer side of the interface and the inner core side respectively. The interface conditions imposed herein thereby eliminate some certain surface instabilities which is a possible when the interfacial adhesion is weak. Alternative treatments that seek to account for the weak or otherwise imperfect bonding conditions that lead to various interface phenomena (e.g., relative sliding displacement or normal direction delamination) which are beyond the scope of this work. For perfect bonding, the traction vector is continuous at the interface T|+ ¯n = T|− ¯n where ¯n is the unit normal to the interface boundary in the deformed configuration 𝜕Ω𝐼 𝑥. In the reference configuration prior to growth and deformation, the normal to the interface is N = E𝑌 . The normal vectors n and N are related by Nanson’s formula n𝑑𝑎 = 𝐽F−𝑇 N𝑑𝐴. The corresponding first Piola- Kirchhoff stress tensor P, especially for later developments, must thus satisfy P N 𝑑𝐴 = T n 𝑑𝑎 for any area element, where N is the unit vector normal to the area element 𝑑𝐴 in the reference configuration. The stress measure P is thus expressed by P = 𝐽TF−𝑇 . (5.2.8) Thus, the interface condition can be expressed as P|+ ¯N = P|− ¯N where ¯N is the unit normal to the interface boundary 𝜕Ω𝐼 𝑋 in the reference configuration. Traction continuity is a natural interface and so does not require explicit enforcement in the finite element procedure. 103 5.2.5 Constitutive function for the energy In this specific problem at hand, a compressible neo-Hookean constitutive model is considered for both layers. The compressibility is accounted for by a direct additive contribution of the pure volumetric part 𝑊VOL(𝐽𝑒) into the deviatoric part 𝑊DEV(𝐼 𝑒 1 ). The volumetric part can be expressed using various empirical laws that lead to different stored energy density models [96]. The attention is confined to the following form of compressible neo-Hookean type stored energy given by (cid:40) 𝑊𝑐 = 𝜇𝑐 2 (𝐼 𝑒 1 − 2 − 2 ln 𝐽𝑒) + (cid:0) 𝜅𝑐 2 − 𝑊 = 𝑊𝑠 = 𝜇𝑠 2 (𝐼 𝑒 1 − 2 − 2 ln 𝐽𝑒) + (cid:0) 𝜅𝑠 2 − 𝜇𝑐 3 𝜇𝑠 3 (cid:1) (ln 𝐽𝑒)2, on Ω𝑐 (cid:1) (ln 𝐽𝑒)2, on Ω𝑠 (5.2.9) where 𝜇𝑐 and 𝜇𝑠 are the shear modulus while 𝜅𝑐 and 𝜅𝑠 represents the bulk modulus of the cortex and subcortex, respectively. In making the correlation to the small strain linear theory, 𝜅 and 𝜇 correlate with Poisson ratio 𝜈 = (3𝜅 − 2𝜇)/(2𝜇 + 6𝜅). For the purposes of the present work, 𝜅 and 𝜇 for both the cortex and subcortex are chosen so as to give a linear theory Poisson ratio of 0.45. This makes (cid:40) 𝑊𝑐 = 𝜇𝑐 2 (𝐼 𝑒 1 − 2 − 2 ln 𝐽𝑒 + 9(ln 𝐽𝑒)2), on Ω𝑐 𝑊 = 𝑊𝑠 = 𝜇𝑠 2 (𝐼 𝑒 1 − 2 − 2 ln 𝐽𝑒 + 9(ln 𝐽𝑒)2), on Ω𝑠 and the Cauchy stress ( 5.1.2 ) thereby becomes T = 𝜇 𝐽𝑒 (cid:0)B𝑒 + (9(ln 𝐽𝑒) − 1)I(cid:1), (5.2.10) (5.2.11) where 𝜇 is the shear modulus of the material i.e., either 𝜇𝑐 or 𝜇𝑠 depending on the location. 5.2.6 Computational implementation The height of subcortex 𝐻𝑠 is ten times bigger than the height of cortex layer 𝐻𝑐 to minimize the effect of the inner core thickness on the surface pattern, i.e., 𝐻𝑠 = 10𝐻𝑐. The domain was then discretized using structured triangular elements using the open-source mesh tool GMSH. The mesh was refined in regions where gradients were expected to be relatively stronger, such as the outer layer and the interface between the outer layer and the inner core. 104 To avoid self-contact at the sulcal gap, we employ a half-wavelength-length domain approach with the selection of suitable displacement bounds to simulate only one-half of the gyrus-sulcus formation. The initial wavelength of the wrinkle 𝜆 is computed using (1.2.1). Note that the film and substrate stiffness in (1.2.1) correspond to the cortex 𝜇𝑐 and the subcortex 𝜇𝑠 shear moduli in the presented context. The length of the domain thereby becomes 𝐿 = 𝜆/2. With the use of suitable upper and lower bounds for each directions, (5.2.4) turns into a bound-constrained minimization problem. The energetically favorable displacement field for each equilibrium was computed using Trust- Region Optimization (TRON) algorithm [82] available in PETSc. The numerical simulation settings were implemented using the open-source scientific computing platform FEniCS [5]. 5.3 Qualitative morphological change and quantitative triggering by the increasing mis- match strain In the presented model, the stiffness ratio 𝜇𝑐/𝜇𝑠 is the key determinant of how the surface patterns emerge following to a certain sequence of transformations as the mismatch strain increases. For a validation purposes, we direct our attention to a specific scenario wherein the stiffness contrasts 𝜇c/𝜇s are 3.64 and 67.24. Figure 5.2 demonstrates how the equilibrium configuration of two separate patterns of morphological development evolves as the mismatch strain increases. At relatively modest values of the mismatch strain 𝜀𝑀, the cortical structure maintains its flat morphology, continuing to expand perpendicularly to the original outer surface. As mismatch strain reaches to a critical value of 𝜀𝑀, the initially flat cortex undergoes a transformation into a wrinkled state. The critical value of 𝜀𝑀 at which this transformation occurs is found to be 0.24 and 0.11 for the respective cases. For the scenario with 𝜇c/𝜇s = 3.64, the wrinkled state starts folding inward at approximately 𝜀 = 0.28 which establish self-contact, while for the case with 𝜇c/𝜇s = 67.24, the transition to period-doubling takes place around 𝜀 = 0.31. As the strain mismatch further increases, amplitude and wavelength of the surface pattern alters. The evolution of the folding patterns and transformations presented above is consistent with the experimental results presented by [136] (cf. Figure 3). It is essential to highlight that the model presented herein is not limited to the surface 105 Figure 5.2 Evolution of the equilibrium configuration for 𝜇c/𝜇s = 3.64 (top) and 𝜇c/𝜇s = 67.24 (bottom) from flat to folding (top panel) and flat to period-doubling (bottom panel) as the mismatch strain 𝜀M increases. patterns demonstrated. Alternative surface patterns, such as creases, can also be formed through the utilization of this model, depending on the stiffness contrast ratio. Specifically, the occurrence of creasing instability, observed when 𝜇c/𝜇s is relatively lower, necessitates specialized treatment, and is notably sensitive to artificial defects. Yet, this instability is beyond this study. Additionally, when considering scenarios involving less-than-perfect bonding, the inclusion of adhesion energy within the model description may be required. This inclusion may result in the emergence of additional types of surface patterns, such as delamination buckling. However, it is noteworthy that this particular instability does not align with the anatomical characteristics of the actual brain. We confine our attention to the folding instability that resembles sulcus-gyrus formation in the context of brain anatomy. Left panel of figure 5.3 presents a comprehensive depiction of the morphological evolution of sulcus-gyrus formation, given in the top row of Fig. 5.2, behave depending on the mismatch strain. In the early stage, the cortex initially runs in the direction perpendicular to the pial surface and remains flat until 𝜀𝑀 ≈ 0.24 (inset a). Afterward, the domain transforms into a wrinkling type instability until the initial fold starts at 𝜀𝑀 ≈ 0.28 where the self-contact takes place (inset b). Throughout the folding phase, the sulcal floor deepens by moving perpendicularly to the initial cortical surface, whereas the gyrus increases in height. The gyrus continues to develop asymmetrically while bringing the adjacent gyral crowns closer together. The ventral wall of the gyrus slopes steeply towards the sulcal floor, whereas the dorsal wall of the gyrus shelves gently into the intrasulcal gap. As growth progresses, a shorter-wavelength secondary fold 106 Folding (Wrinkle with cusped end)WrinkleSecondary Fold=0.24=0.33=0.41=0.44=0.48FlatεMεMεMεMεMWrinklePeriod-DoublingεM=0.24=0.31=0.11=0.48=0.58εMεMεMεM Figure 5.3 Evolution of the morphological features of a sulcus-gyrus formation as the mismatched strain increases. Scalings to evaluate morphological features of the final folding formation [125]: (− · · −) 𝜆/𝐻𝑐, (− · −) 𝑑/𝐻𝑐, and (—) ℎ𝑔/ℎ𝑠. (- - -) where 𝜆: the distance between two adjacent gyral crown, 𝐻𝑐: initial thickness of the cortex, 𝑑: distance between the top of the gyral crown to the deepest point of the sulcus gap, ℎ𝑔: thickness of sulcus, as illustrated in the inset. Strain mismatches: 𝜀M = 0.24 (flat to wrinkle), 𝜀M = 0.28 (wrinkle to fold) and 𝜀M = 0.44 (shorter-wavelength secondary folding)(Left panel). Change of self-contact length as the mismatched strain increases. (−) the self-contact length showed in the insets. (Right panel). thickness of gyrus, and ℎ𝑠: forms around 𝜀𝑀 = 0.44 (inset c), increasing the wavelength of the folding nearly twice. The dorsal part of the gyrus continues to grow outwards and becomes roughly symmetrical. The cortex thickness also varies, resulting in the gyrus becoming thicker than the cortical plate in the sulcus during development. These findings are consistent with (i) the numerical model results and the measurements obtained from the porcupine, cat, and human brains for the same scaled dimensions displayed in Fig. 2 of [125], (ii) structural changes observed in postnatal ferret brain experimentally [119]. The right panel of Figure 5.3 depicts the evolving characteristics of the self-contact area in response to increasing strain mismatch. This self-contact region is analogous to the sulcal gap within the context of cerebral anatomy. Specific troughs within the wrinkles undergo an inward folding process, resulting in the establishment of self-contacts. This folding occurs after transitioning from the wrinkled configuration, typically at a critical strain value of approximately approximately 𝜀𝑀 ≈ 0.28. Consequently, two points on the outer surface of the cortex, initially separated by a finite distance, first make contact with each other. As the strain mismatch is further elevated, the 107 h(c)(a)dghsFlat WrinkleFoldingSecondary Fold(b)FlatWrinkleFoldingSecondaryFoldd/Hchg/hsλ/Hcλ self-contact, which was originally at a point, becomes a line of self contact. The length of this self contact line increases with 𝜀𝑚 until secondary folding occurs. For the process depicted in left panel of Fig. (5.3) this secondary folding initiates at a critical strain of approximately 𝜀𝑚 = 0.44. Just prior to secondary folding, the original contact line is at its maximum length. Upon secondary contact, the original contact line abruptly decreases in length as certain points previously in contact now separate. This is evident in right panel of Fig. (5.3) where the original maximum contact length of 2.11 is suddenly reduced to a length of 1.24. This secondary folding process gives rise to a surface pattern characterized by a shorter wavelength and reduced amplitude. 5.4 Residual Fields in the Sulcus-Gyrus Formation 5.4.1 Residual strain fields In conjunction with the morphological alterations, slow and substantial deformations arising from growth and gyrification exert a significant influence on the internal mechanical state of the tissue, thereby playing a crucial role in comprehending the homeostasis conditions within the natural state of the tissue. Therefore, the attention is now directed to spatiotemporal evolution of residual stress and strain fields within a sulcus-gyrus formation. To this end, the Green-Lagrange strain measure is employed to assess the large deformations emerging during the flat-to-folding transformation presented above, i.e., 𝜺 = 1 2 (C − I), where C is the right Cauchy-Green tensor C = F𝑇 F. It should be noted that 𝜺 is related to the total deformation meaning that 𝜺 includes the growth and the effect of elastic deformation. The elastic deformation strain is expressed by 𝜺𝑒 = 1 2 (C𝑒 − I), where C𝑒 is given in (1.3.5). Because of the absence of the growth strain in the subcortex, the total strain only includes the elastic effects i.e., 𝜺 = 𝜺𝑒. However, the total strain 𝜺 is the combination of both the elastic strain and growth strain together in the cortex. The tangential growth in cortex is constrained by the non-growing subcortex. During the transition from a flat to a wrinkled state occurring at approximately 𝜀𝑀 = 0.24, the tensile elastic strain manifests in the cortex and certain regions of the subcortex that corresponds to the crest of the wrinkled state in the direction of cortex folding. As the gyral crown experiences an increase in height, the tensile elastic strain becomes notably localized within the subcortex region 108 of the gyrus, specifically emphasizing its presence at the cortex-subcortex interface. With the initiation of folding (0.27 < 𝜀𝑀 < 0.45), the accumulation of compressive (negative) elastic strain commences within the regions that eventually constitute the sulcal floor as the folding progresses. At the mature sulcus-gyrus formation, compressive elastic strains become evident within the depths of the sulcus and the subcortical areas of the sulcus, particularly along the gray-white matter interface. Simultaneously, the substantial tensile residual strain becomes apparent throughout the entire sulcus and gyrus, with localized concentrations observed at the sulcal floor and the pial surface of the gyral crown, propagating in the growth direction. Notably, the magnitude of tensile strain in these regions approximates tenfold the maximum recorded compressive strain observed within the subcortical domain of the cortex. These results are displayed in figure 5.4. Figure 5.4 Evolution of the residual strain fields in 𝜀𝑥𝑥 (column 1, 2) and 𝜀𝑦𝑦 (column 3, 4) at four different strain mismatch states: 0.24, 0.33, 0.41, 0.41, and 0.51% (from top to bottom). The strain fields are isolated based on the sign: compressive (column 1, 3) and tensile (column 2, 4). The strain has opposite nature in non-colored areas. 109 1.8-0.45-0.353.02.01.0-0.4-0.20.00.00.51.01.50.00.0εyyεyyεxxεxx3.7-0.3-0.2-0.1 5.4.2 Residual stress fields in sulcus-gyrus formation: the first PK stress Figure (5.5 ) depicts the contours of the isolated residual stress field during the sulcus-gyrus formation. To assess the stress field, the first Piola-Kirchhoff stress measure P given in (3.8) is employed to map the spatial distribution of the stress in the x- (tangential to the reference cortical surface) and y-directions (normal to the reference cortical surface) as illustrated left and right two panels respectively. The results reveal the existence of growth-induced residual stress fields that are primarily localized in the cortex and the regions of the subcortex particularly close to the interface. Both the cortex and subcortex exhibit inhomogeneous residual stress fields, which vary both spatially and temporally. Notably, the cortex demonstrates a compressional residual stress field in the x-direction maintaining this consistent behavior throughout the transformation process. The compressional 𝑃𝑥𝑥 is then localized at the innermost part of the mature gyrus. It should be noted that only the pial (outermost) surface of the gyrus exhibits tensile stresses in this direction. The results further reveal a notable transition from compression to tensile residual stresses at the gray and white matter interface, occurring in both the gyral crown and the sulcal floor. Along with the pial surface of the gyrus, two specific regions within the subcortex exhibit tensile residual stress fields: the inner part of the gyrus and the sulcal floor. In the y-direction, as shown in the right two panels, there is a predominance of compressional residual stress at the sulcal floor, concentrated around the deepest point of the sulcal gap. Conversely, the sulcal walls consist of the tensile residual stress field along the same y-direction, which exists around the self-contact region of the folding during throughout the transformation. The subcortex region of gyrus is the other regions experiencing the tensile stress in y-direction. The reference configuration normal to the outer surface and to the interface is N = E𝑌 . The associated reference configuration traction vector is PE𝑌 = 𝑃𝑦𝑥E𝑋 + 𝑃𝑦𝑦E𝑌 . The free surface boundary condition discussed in Section 3.3 thus gives that 𝑃𝑦𝑥 and 𝑃𝑦𝑦 must vanish at the top surface. Due to the interface traction continuity i.e., P|+E𝑌 = P|−E𝑌 , the components 𝑃𝑦𝑦 and 𝑃𝑦𝑥 are continuous at the interface. These observations are consistent with figure 5.5. 110 Figure 5.5 Tensile (first and third column) and compressional (second and fourth column) residual stress fields in x-(left two panels) and y-(right two panels) at four different strain mismatch states: 0.24, 0.33, 0.41, and 0.51 (from top to bottom). The stress has opposite nature in non-colored areas. 5.4.3 Residual stress fields in sulcus-gyrus formation: Cauchy Stress Figure (5.6) depicts the contours of the isolated residual stress field during the sulcus-gyrus formation. To assess the stress field, the first Cauchy stress measure 𝝈 given in (5.2.11) is employed to map the spatial distribution of the stress in the x- (tangential to the reference cortical surface) and y-directions (normal to the reference cortical surface) as illustrated left and right two panels respectively. The results reveal the existence of growth-induced residual stress fields that are primarily localized in the sulcus floor, deep gyrus and the regions of the subcortex close to the interface. Both the cortex and subcortex exhibit inhomogeneous residual stress fields, which vary both spatially and temporally. Notably, the cortex demonstrates a compressional residual stress field in the x-and y-directions maintaining this consistent behavior throughout the transformation process. The compressional 𝜎𝑥𝑥 is then localized at the innermost part of the mature gyrus and the deepest point of sulcus. It should be noted that only the pial surface of the gyrus exhibits tensile stresses in this direction. The results further reveal a notable transition from compression to tensile residual stresses at the gray and white matter interface, occurring in both the gyral crown and the 111 01230-10-28-200-20-40-50PxxPxxPyyPyy1202.5μsμsμsμs sulcal fundus. In addition to the pial surface of the gyrus, the gyral white matter and the subcortex part of sulcal fundus are two other specific regions existing tensile stress field. In the y-direction, as shown in the right two panels, there is a predominance of compressional residual stress at the sulcal fundus, concentrated around the deepest point of the sulcus. The subcortex region of the gyrus is the region experiencing the tensile residual stress in the y-direction. Figure 5.6 Tensile (first and third column) and compressional (second and fourth column) residual stress fields in x-(left two panels) and y-(right two panels) at four different strain mismatch states: 0.24, 0.33, 0.41, and 0.51 (from top to bottom). The stress has opposite nature in non-colored areas. We make use of octahedral shear stress theory to evaluate how distortional strain energy changes within the tissue during transformations. To this end, we utilized the von Mises stress which is related to the second invariant of the deviatoric stress s, i.e., 𝜎v = ( 3 2s : s)1/2. Here, the deviatoric 𝛿𝑖 𝑗 𝜎𝑘 𝑘 , where 𝝈 denotes the Cauchy stress tensor. As depicted stress s is expressed by 𝑠𝑖 𝑗 = 𝜎𝑖 𝑗 − 1 2 in Figure (5.7), the von Mises stress is prominently observed in the cortex during development. The von Mises stress accumulates specifically around the innermost part of the future gyrus and around the sulcal walls as the mismatch strain elevates. In the mature gyrus-sulcus formation 𝜀𝑀 = 0.51, the maximum von Mises stress is observed in the innermost region of gyral gray matter particularly 112 0-10-20-300246.60-10-20-3002.254.5σxxμsσxxμsσyyμsσyyμs close to the cortex-subcortex interface. Likewise, it is noteworthy that the sulcal floor represents another spatial domain where von Mises stress attains a notably elevated state. These regions exhibit a susceptibility to undergoing permanent deformations. In accordance with the deformation theory of plasticity, such deformations could give rise to a nonlinear constitutive relation and a softened tangent modulus in these specific regions with the onset of the permanent deformations process. It is experimentally shown that brain tissue exhibit Mullins effect and hysteresis behavior [42]. Figure 5.7 The von Mises stress field during the gyrus-sulcus formation at different strain mismatch levels. It should be noted that every residually-stressed body has a non-zero strain energy density field in its deformed configurations. Hence, the strain energy density function can be utilized as a metric to quantify the residual stress field [75]. Figure (5.8) displays how the strain energy density function evolves within the material during the gyrus-sulcus formation. As the mismatch strain increases, the strain energy accumulates specifically around the innermost part of the future gyrus and around sulcal fundus. Figure 5.8 Evolution of strain energy density field during the gyrus-sulcus formation for given mismatch strains 𝜀𝑀. 113 50101520εM=0.24εM=0.33εM=0.44εM=0.51σVμs2045εM=0.24εM=0.33εM=0.44εM=0.51Wμs 5.4.4 Residual stress field distributions in the sulcal fundus, gyral crown and sulcal wall As evident from the top panel of figure (5.9), the residual strain fields progressively increase along the cortex until they reach the interface. The maximum strains are measured at the interface, where 𝜀𝑥𝑥 attains a magnitude nearly 6 times higher compared to 𝜀𝑦𝑦. While 𝜀𝑥𝑥 is elongational, 𝜀𝑦𝑦 has a contractile nature along the cortex and interface. Subsequently, the both strains diminish along both axes as they extend into the bulk of the tissue. In contrast, the maximum stress is observed at the tip of the sulcus for both axes, gradually relaxing in an asymptotic manner as it approaches the interface. At the interface, the normal Cauchy stress component 𝜎𝑥𝑥 and 𝜎𝑦𝑦 exhibit a slope discontinuity. While 𝜎𝑥𝑥 undergoes a transition to the tensile field and then continues to relax asymptotically towards zero stress, 𝜎𝑦𝑦 shows a similar behavior by remaining compressive. In the bottom panel of figure (5.9), the highest values of 𝜀𝑥𝑥 and 𝜀𝑦𝑦 are recorded at the pial surface and the gray-white matter interface, respectively. However, the discontinuity in 𝜀𝑦𝑦 is not as strong as that observed in the mature sulcus, and it gradually reduces in the subcortex, displaying a relatively linear trend. Similarly, 𝜀𝑥𝑥 diminishes along the cortex with a relatively linear trend and follows a steady behavior in the subcortex. The results also suggest that the residual strains in the mature sulcus exhibit the same sign in each layer, whereas the residual strains encompass both positive and negative values in the mature gyrus. Regarding the normal stress components 𝜎𝑥𝑥 and 𝜎𝑦𝑦 progressively increases until it reaches the maximum magnitude of compressive stress at the interface. Notably, 𝜎𝑥𝑥 undergoes a relatively strong discontinuity compared to 𝜎𝑦𝑦, subsequently reversing and decaying towards zero stress. Conversely, the gyral white matter contains a tensile residual stress in the y-direction 𝜎𝑦𝑦 with moderate value. In conjunction with the analysis of stress fields within the sulcal floor and the gyral crown, we also extend our attention to the stress distribution in the sulcal wall. The formation of the sulcal wall arises from the transition from a wrinkled state to a folded state, a transformation occurring at a critical strain magnitude 𝜀𝑀 ≈ 0.27. This transition emerges as the convergence of distinct points on the outer cortical surface, originally separated by a finite spatial interval, thereby forming a contact line corresponds to the sulcal gap. As the transformation progresses, depending on the 114 Figure 5.9 Green-Lagrange strain (left column) and Cauchy stress (right column) distributions in the 𝑥- (dashed blue) and 𝑦- (red solid) directions are shown along the probe line for the sulcus (top panel) and gyrus (bottom panel). For the sulcus, the probe line starts from the deepest point of the intrasulcal gap and extends through the gray matter (GM) and white matter (WM) respectively. For the gyrus evaluation, the probe line starts from the pial surface of the gyral crown and extend towards gyral white matter. The location of the probe line is also shown in the associated insets. The dashed black line represents the gray-white matter interface location. level of strain mismatch, the contact line undergoes dynamic alterations, as previously illustrated in figure 5.3. The residual stress field in a developing sulcus-gyrus formation was also previously demonstrated via the contours of figure 5.6. It indicates the presence of a compressive residual stress field in the x- and y- direction. Figure (5.10) presents the spatial distribution of residual stress along the defined probe lines given in corresponding insets. The placement of probe lines is strategically aligned along the contact region of the sulcal wall. Namely, the probe line in the left panel of Figure (5.10) starts at the point of contact initiation and extends towards to the deepest point of the sulcus. Simultaneously, the probe line in the right panel commences from the midpoint of the sulcal wall, extending horizontally into 115 GWMyx001.33GMWMyx02.285.0 the gyral white matter. The results indicate a compressional stress field for both normal components of Cauchy stress 𝜎𝑥𝑥 and 𝜎𝑦𝑦, exhibiting an increment along the contact line and concentrating at the deepest point of sulcus. Likewise, a compressional residual stress field is observed along the probe line depicted in the right panel inset. Notably, upon reaching the interface, 𝜎𝑥𝑥 exhibits a discontinuity and rapidly drops to markedly low magnitude, while 𝜎𝑦𝑦 turns into a tensile stress field and continues to exist with a relatively moderate magnitude. Figure 5.10 Cauchy stress distributions along the probe line illustrated in the associated insets. (· · ·) and (−−)the deepest of sulcus, and the location of the gray-white matter interface, respectively. Panel b of Figure 5.3 showed how folding gave rise to self contact along the sulcal wall, the length of which increased with the mismatch strain 𝜀𝑀. This contact length was abruptly reduced by the advent of secondary folding, after which the contact length again increased with 𝜀𝑀 . The two contact sides push against each other with traction 𝝈E𝑋 = 𝜎𝑥𝑥E𝑋 + 𝜎𝑥𝑦E𝑌 . Figure (5.11) shows the evolution of the normal contact stress 𝜎𝑥𝑥 with mismatch strain 𝜀𝑀 during the first folding regime 0.27 ≤ 𝜀𝑀 ≤ 0.45 (left panel) and the subsequent secondary folding regime 0.45 ≤ 𝜀𝑀 ≤ 0.51 (right panel). It should be noted that the maximum point of the contact length corresponds to the deepest point of the sulcus in both panels. In both regimes, the normal contact stress is compressive along the self-contact line with monotonically increasing magnitude along the contact line. Regardless of the mismatch strain level, the deepest point of the sulcus displays a highly localized maximum contact stress that consistently displays values close to -30𝜇𝑠. 116 GMWMyx01.64Deepest Point of sulcusGMWMyx0230 Figure 5.11 Distribution of normal contact stress along the line of self-contact in the primary fold as determined by the value of mismatch strain. The left panel depicts the range 0.27 ≤ 𝜀𝑀 ≤ 0.45 for primary folding, and the right panel depicts the range 0.45 ≤ 𝜀𝑀 ≤ 0.51 in which a secondary fold is developing between the primary folds. 5.4.5 Analysis of principal stress and directions The normal components of stress tensors do not inherently align with the resultant morphology of the gyrus-sulcus formation. To address this, the investigation is extend into the examination of the eigenvalues of the Cauchy stress tensor and the eigenvector associated with these eigenvalues, corresponding to the principal stresses and principal directions. These are given by the two eigenvalues and two eigenvectors associated with the in-plane part of the stress tensor at hand. The z-direction serves as the third principal direction. The stress component 𝜎𝑧𝑧 acts as a reactive wall stress in the presented plane strain model, thereby inducing a compressive stress field along this direction. Due to the compressive nature, the eigenvalue associated with the z-direction would correspond to either intermediate or minimum principal stress in the full 3-D treatment. Locations where the maximum principal stress from the 2-D treatment is negative correspond to places where there is no tensile aspect to the stress field. Such locations are relatively resistant to crack formation within the tissue. Our special interest is in locations where the maximum principal stress from the 2-D treatment is positive, since the resulting tensile stresses can promote crack formation and other possible failure behavior of the tissue. For this reason, only the evolution of positive (tensile) maximum principal stress is depicted in the contours of figure 5.12. During the wrinkling phase shown at 𝜀𝑀 = 0.24, the maximum principal stress in tension primarily localizes 117 0.450.320.410.270.370.450.460.490.480.51 at the lateral side of the interface, with the exception of the future sulcal floor. Subsequently, as the folding is initiated at 𝜀𝑀 = 0.27 and deepens through the subcortex 0.27 ≤ 𝜀𝑀 ≤ 0.45, the principal stress fields begin to accumulate prominently along the sulcal floor and certain regions of the cortex that will eventually transform into gyrus. Upon attaining the mature gyrus-sulcus formation, it becomes evident that the maximum principal stresses in tension concentrate predominantly in two specific regions: (i) along the gray-white matter interface along the sulcal floor and (ii) at the tip of the gyrus within the subcortex. The reduced dimension of the stress tensor yields two analogous eigenvectors associated with the maximum and minimum eigenvalues thereof. These eigenvectors represent the principal directions along which the maximum and minimum principal stresses act within the gyrus-sulcus formation. However, our focus is confined to only the eigenvector associated with the maximum (positive) eigenvalue since any crack opening is driven by the maximum tensile stress acting along the associated direction. In the context of biomechanics, the term “fissure" more aptly captures the essence of a crack within a soft tissue material. The direction of any such hypothetical fissure is oriented at a 90◦ from the direction of the maximum principal stress. In conjunction with the maximum (positive) principal stress contours, the direction of the possible hypothetical fissure opening is also presented in figure 5.12 by the black lines. The length of the lines is proportional to the magnitude of the maximum principal stress. The fissure trajectories exhibit a strong perpendicularity to the interface where the maximum tensile principal stress field is observed. As the analysis progresses towards regions with maximum tensile stresses with lower magnitudes (depicted by blue color), the orientations of the fissure apertures evince a gradual convergence towards parallelism with the fold structure. 118 Figure 5.12 Evolution of the maximum principal stresses for only tensile (positive) fields and corresponding fissure pattern (dashed black lines) within the gyrus-sulcus formation. The principal stress has a reverse sign in non-colored area. 5.4.6 Effect of skull constraint on the morphology and residual fields Previous experimental work shows that the constraint generated by the surrounding intracranial structures is not necessary to trigger the cortical folding [8]. Still, the presence of a skull can influence the morphology of the cortex, particularly by flattening the gyral crown. In this section, the effect of skull constraint on the sulcus-gyrus formations and on the mechanical field in the resultant folding pattern is investigated. In the computational model, the confinement is implemented through the upper bound of the displacement field in the y-direction, i.e., 0 ≤ 𝑢𝑦 ≤ 𝐻 + 𝐷 where 𝐻 is the initial height of the domain, and 𝐷 is the distance between the outermost surface of the cortex and the skull. The presence of a skull here acts as a boundary condition that constrains the 119 0.60.40.20.80.06.04.02.06.70.01.00.51.30.04.02.05.40.01.51.00.52.20.03.02.01.04.20.0εM=0.24εM=0.28εM=0.33εM=0.41εM=0.44εM=0.51σμsσμsσμsσμsσμsσμs displacement in a direction normal to the boundary. The level of confinement is arranged by means of the parameter 𝐷 while the other growth and boundary conditions presented above remain the same. Figure 5.13 shows how the skull constraint affects the resultant sulcus-gyrus formation and corresponding stress contour at the strain mismatch 𝜀𝑀 = 0.51. In the computational model at hand, the parameter 𝐷 was assumed as a sufficiently high value to exclude the effect of skull constraint. This case is presented in the bottom row of Figure 5.13. Under this condition, the gyral crowns exhibit a relatively curved outer profile, while the tips of the subplate tend to be relatively cusped. The displacement of the outermost point of the gyrus in the y-direction 𝑢𝑦 was obtained as 2.56 mm at 𝜀𝑀 = 0.51. Hence, all 𝐷 values less than the maximum displacement 𝐷 < 2.56 will automatically lead to constraints in this direction under the same conditions. To test this effect, the parameter 𝐷 is assigned as 2 mm to generate a constraint acting like a skull. As shown in the leftmost panels within the figure, the results reveal that the skull constraint results in flattened regions within the gyrus, notably at two specific loci: the pial surface and the point of junction between the subplate and the white matter. This observation aligns well with prior clinical investigations and corroborates the findings reported in Tallinen et al. (2014). It is noteworthy that the presence of the skull constraint does not exert any influence on the temporal evolution of the gyrus-sulcus formation such as critical strain mismatch. Figure 5.13 Effect of skull constraint on the sulcus-gyrus formations and the normal residual stress fields for 𝐷 = 2 mm (top panel), and 𝐷 = 4 mm (bottom panel) at the strain mismatch 𝜀 = 0.51. 120 0-10-20-300246.60-10-20-3002.254.5σxxμsσxxμsσyyμsσyyμs Furthermore, the results reveal that the skull constraint has the potential to affect the residual stress field within the formation, as seen in the stress contour given in the figure. Overall the distribution of the stress fields shows a certain similarity. Still, the skull constraint leads to notable changes in residual stress fields specifically at the interface of the cortex and subcortex in the gyrus. In the unconstrained model, the stress fields, and in particular, 𝜎𝑥𝑥 exhibiting compression, and 𝜎𝑦𝑦 displaying tension, tend to localize prominently at the cusped tip of the inner gyrus. However, more uniform stress distributions are observed at these locations in the models where the skull constraint is considered. In addition to the stress contours, Figure 5.14 shows the effect of skull constraint on the distribution of residual stress fields along the probe line. The probe line begins from the pial surface of the gyrus and extends through the bulk of the tissue. While the qualitative characteristics of the stress fields, particularly 𝜎𝑥𝑥, remain consistent, there is a notable reduction in the magnitude of compressive stress at the interface in cases where the skull constraint is imposed. Furthermore, in the unconstrained model, the stress field 𝜎𝑦𝑦 exclusively exhibits a single sign (positive), whereas the cortex in the constrained model includes both positive and negative signs. The observed discontinuity at the interface becomes less pronounced for both normal components of the residual stress field. It is observed that the skull constraint has minimal effect on the qualitative behavior of the stress field along the sulcal floor. The distribution of residual stress field along the sulci is not therefore represented intentionally. 121 Figure 5.14 Distribution of the normal components of residual stress field 𝜎𝑥𝑥(blue) and 𝜎𝑦𝑦(red) for constrained (dashed, C) and unconstrained (solid, UC) models.The dotted and dashed lines show the location of the interface for C and UC models respectively. 122 UCUCCC 5.5 Discussion 5.5.1 Observations Previous experimental works [143, 144] have demonstrated the qualitative attributes charac- terizing residual stress fields within mouse and evolving ferret brains. The experimental protocol focused on unveiling the stress distribution of the coronal section of brain tissue by executing radial and tangential incisions with various depths. These incisions engendered stress relaxation orthogonal to the cutting trajectory. In the mouse experiments, the observations are confined to only tangential residual stress distributions due to the direction of incisions. The results indicate that gray matter is in compression whereas the white matter has a tensile residual stress field. Fur- thermore, it is important to note that the mouse brain has a lissencephalic structure, characterized by a lack of convoluted contours on its outer surface, unlike the mammalian brain. Hence, the growth of the mouse brain resembles the early phase of the computational model at hand where the layers expand solely prior to folding i.e., 𝜀𝑀 ≤ 0.24. The dissection experiment performed on the developing ferret brain emerges as a more appro- priate endeavor to understand how the growth and gyrification shape cerebral morphology, and the internal mechanical state of gyral and sulcal formations. It should be noted that the cerebral architecture of the ferret is characterized by a relatively modest degree of folding. The qualitative attributes characterizing residual stress fields within a fully matured ferret brain are summarized in Table 5.1. This table demonstrates the comparative analysis between the outcomes acquired from the presented model and the experimental observation. The consistency between the residual fields intrinsic to the gyrus-sulcus formations and the empirical observations is evident. The sole point of divergence is observed in the deep white matter of the gyrus where the model estimates only compression whereas the experiment indicates a tensile residual stress field. For the same region, the localized tensile residual stress at the gyrus interface predicted by the model have not been observed in the fully mature ferret brain. Further investigation is imperative to examine the dynamics of post-instability regimes, particularly as manifested in extensively convoluted brains such as the human brain. 123 Table 5.1 Comparison of the result of dissection experiment conducted on fully mature ferret brain and estimations of the computational model at hand. The experimental observations have been presented in [144] (cf. Figure 3 c1-c6). G: Gyrus, S:Sulcus, GM: Gray Matter, WM: White Matter, XX: Tangential to the folding direction, YY: Normal to the folding direction. Region Tissue Direction Dissection Experiment Computational Model G G S S G G S S GM WM GM WM GM WM GM WM XX XX XX XX YY YY YY YY Tension → Compression Tension → Compression Compression → Tension Compression Tension Compression Tension → Compression No Data Compression Compression Compression Tension Compression Tension Compression Compression The stress states predicted by the computational model are consistent with previous experimental observations. Additionally, the evolution of equilibrium states of the computational domain is parallel to the gyrus and sulcus formation experimentally observed in [119]. Hence, it is fair to assume that the analysis and simulations presented here support the hypothesis that differential tangential growth of the cortex. The results also show that the growth and gyrification lead to the nonhomogeneous residual fields within the gyrus and sulcus. The sulcal floor and gyral crowns are the regions mostly deformed regions in the gyrus-sulcus formation. Considering the effect of residual stress in the realm of solid mechanics, these mechanical intrinsic factors can help to explain various phenomena regarding heterogeneity in material properties, interface conditions, and injury kinematics. To the author’s knowledge, this is the first detailed investigation addressing growth- induced residual stress (strain) fields within the gyrus-sulcus formation. Previous computational studies have various motivations to investigate various aspects of cerebral growth and folding except for the growth-induced mechanical intrinsic factors. Previous work [125] demonstrated the residual stress field of a single sulcus for validation purposes. The computational work of [9] presented stress distributions in a growing bilayer system; however, the computational predictions are limited only to pre-folding phases, mostly wrinkling. 124 5.5.2 Limitations of the model The computational model provides an intuitive understanding of how the residual stress and strain fields evolve in a developing gyrus-sulcus formation, despite the relatively low level of elaboration of the computational model. Yet, there are multiple limitations to this model. In this work, a specific stiffness ratio 𝜇𝑐/𝜇𝑠 = 3.64 is utilized. It is challenging to characterize the mechanical properties of a developing brain and the substantial changes that brain tissue undergoes at the early stages of the brain. Because of a lack of experimental data, the stiffness modulus contrast at hand could not be fully representative. Recent experimental works suggest that the gray and white matter in the mature brain is equally stiff, i.e., 𝜇𝑐/𝜇𝑠 ≈ 1. However, our current understanding of the material properties of gray and white matter is significantly limited in the early stages of brain development. Given the spatiotemporal alterations that occur in the microstructure of brain tissue during development, it is reasonable to infer that the mechanical properties of gray and white matter shall be distinct from those observed in the mature brain. Specifically, it is proposed that the gray matter may exhibit greater stiffness than white matter due to advancements in microstructural organization leading to increased cross-linking [58]. As the brain undergoes further development, the myelination process and the emergence of astrocytic branches contribute to the increased stiffness of white matter. Therefore, during the initial phases of development, it is conceivable that gray matter may be stiffer than white matter [141]. Secondly, the computational model at hand employed an idealization of a small portion of the brain’s soft tissue. The model does not encompass a detailed anatomical representation of the complete intracranial structures (e.g., skull, CSF, etc.). The presence of the skull and CSF could potentially alter the boundary conditions, exerting discernible pressures upon the evolving structure. While it is established that the skull is not a necessary trigger for cortical folding, as demonstrated by in [8], its presence has a modulatory influence on the resultant morphology, particularly through its propensity to flatten out the curvature of gyral crowns [125]. Furthermore, the role of the skull extends beyond mere shaping the cortex and would have a role in cortical folding dynamics [90]. In our computational model, the heterogeneous configurations inherent in the cerebral anatomy are 125 not include in the model either. Although experiments have shown that the cortical thickness is no longer homogeneous, even at a young age, in humans and even ferrets [105]. An isotropic neo-Hookean constitutive model has been utilized to simulate the response of the tissue to growth and folding. However, the brain tissue exhibits anisotropy, viscosity, permanent deformation, hysteresis, and a highly nonlinear relationship between stress and strain [42]. In the presence of viscoelasticity, the transient nature of the problem, as well as the geometric and material nonlinearities, could lead to multiple stable equilibria, thereby affecting the mechanical response of the material, and ultimate equilibrium state [151]. Furthermore, experiments to characterize the stiffness of the brain are rare, and measured stiffness values vary in a large range. Considering the effect of stiffness contrast on the folding pattern, it is imperative to underscore that the scarcity in tissues’ stiffness quantification, particularly during gestational phases, curtails the precision of parameter selection. 5.5.3 Conclusions In this chapter, an elementary computational model is utilized to examine how mechanical stress builds up in a developing gyrus-sulcus formation. The tangential expansion of the cortex leads to a folding pattern with spatially varying residual stress and strain fields. The resultant folding pattern resembles the characteristic structures of a gyrified brain: sulcus and gyrus. Notably, the computational findings demonstrate qualitative agreement with the results obtained from dissection experiments conducted on the evolving ferret brain. Future studies in this field could explore several intriguing avenues. Firstly, a parametric study can be considered to investigate the sensitivity of the mechanical state emerging within the sulcus- gyrus formation to various parameters. The potential parameter space may consist of: (i) The stiffness modulus contrast 𝜇𝑐/𝜇𝑠. As presented in figure 5.2, the stiffness modulus is one of the key determinants shaping the surface patterns that manifest through a specific sequence of transformations. While the material typically exhibits sinusoidal wrinkling at low levels of mismatch strain, variations in the stiffness modulus ratio, whether higher or lower, give rise to distinct surface instabilities such as folding, period doubling, and period-tripling as the degree 126 of mismatch strain increases, as discussed by Wang [137]. Not only the resultant morphology, recalling (5.2.11), it is obvious that the stiffness of the layers is also effective on the resultant stress field so it needs to be investigated further. In this work, the primary goal is to the examination of the residual fields in a folding instability resembling sulcus-gyrus formation for a specific stiffness ratio 𝜇𝑐/𝜇𝑠 = 3.64. Previous experimen- tal and computational results show how the resultant folding pattern is sensitive to the stiffness ratio. Consequently, it is imperative to comprehensively understand the influence of the stiffness ratio on the emergent residual stress fields within the formation. To conduct such a parametric investigation, it is significant to select the range of stiffness ratios under consideration. For instance, the stiffness modulus less than 1.3 can lead to a crease formation which is different than the wrinkling i.e., 𝜇𝑐/𝜇𝑠 ≤ 1.3. The stiffness modulus ratio should be therefore in the range of 1.3 ≤ 𝜇𝑐/𝜇𝑠 ≤ 5 to investigate the residual fields in folding. (ii) The initial thickness of the cortex 𝐻𝑐 and the distribution of thickness along the subcortex: Considering the equation (1.2.1), it is evident that the thickness of the cortex prior to any growth and gyrification affects the wavelength of the initial sinusoidal wrinkling pattern. This observation aligns with previous computational works that have investigated the influence of initial cortical thickness on the emerging folding pattern [135, 19]. Specifically, a cortex with an initially reduced thickness tends to manifest numerous smaller gyri and shallower sulci, a morphological configura- tion akin to the features observed in polymicrogyria. Along with the thickness of the cortex, how the thickness of the cortex varies along the subcortex is significant and needs to examined further. Namely, the presence of inhomogeneities in cortical thickness could lead to the development of more complex sulcus-gyrus morphologies by acting like a mechanical perturbation [33]. Despite the role of the initial thickness on the morphology, (5.2.11) implies that the initial thickness may not be as much as explicitly effective on the residual stress field. (iii) The growth difference between the layers: In order to emphasize the tangential growth dif- ferences, our consideration is confined solely to the expansion of the cortical surface tangentially, while the white matter does not undergo any concurrent growth. However, the subcortical volume 127 experiences an initial accelerated expansion trajectory, subsequently yielding precedence to the cortical volume as gestational progress ensues, in reality, [7]. It can be assumed that the growth of white matter may be overestimated whereas the growth of gray matter could be overestimated in the presented computational model. Along with that, a spatially uniform and time-independent growth parameter 𝑔 is utilized, which controls the amount of tangential growth, in the computa- tional model. Time-dependent growth functions motivated by the biological changes seen during the morphogenesis offers an interesting research question. Secondly, further work is required to enhance the presented computational model for some aspects. Specifically, augmenting the model to include time-dependent (viscous) effects could furnish a more advanced depiction of the tissue’s growth response (e.g., a visco-hyperelastic constitutive form). There exists a need for further refinement, particularly with respect to the incorporation of the spatial distribution of mechanical and geometrical properties of the tissue, thus transcending the utilization of bulk tissue properties. It is noteworthy that such inhomogeneities in tissue could also potentially affect the surface topology by acting as a perturbation. Additionally, integrating a growth function that accounts for biological changes within a higher dimensional representation would contribute to a more comprehensive simulation. The effect of modifications in the parameter space of this model on growth and gyrification remains an open question and requires further investigation, as well. Lastly, the existence of residual stress fields holds the potential to significantly influence the response of the brain soft tissue and the resulting injury pattern when subjected to an applied load, such as impact. Additional commenting on the possible role of residual stress on brain dynamics is discussed in Chapter 6. 128 CHAPTER 6 CONCLUDING REMARKS, SPECULATIONS, BROADER CONNECTIONS AND FUTURE DIRECTIONS This research explores how a residual stress field can be treated in the framework of hyperelastic theory. The cause of the residual stress field need not be specified during the analysis and, as shown in Chapter 2, a relatively simple nontrivial form for the stress field can yield useful insights. Subsequently, attention is restricted to a residual stress field induced by the differential growth in adjoining tissue. To this end, a numerical analysis was conducted to compute the residual stress fields in the spherical shells for different growth conditions. A finite element model is then built to examine the residual stress fields resulting from an extensive growth period including symmetry- breaking bifurcations. Lastly, the response of single sulci to cavitation-induced deformations in the intrasulcal region was presented, and the potential link between the cavitation deformations and neuropathies/tauopathies such as CTE was discussed. In the second chapter, a new example of a solution of a finite deformation boundary-value prob- lem for a residually-stressed spherical shell subject to pressure-inflation deformation is provided. This study offers obvious avenues to explore further. Following the presented framework, one can consider extending the effect of residual fields analysis for (i) more complicated deformation states, and (ii) different constitutive models rather than the well-known Mooney-Rivlin. On the other hand, it is important to remark that the presented model may not offer a comprehensive representation of the effect of residual stress. One of the reasons is that a simple nontrivial stress field in the radial direction was considered in the study due to the lack of quantitative data at present. The quantitative behavior of the residual stress field needs to be known to develop more elaborate models. The limited data on the state of stress motivates examining the residual stress fields itself arose from a specific reason such as growth. The work therefore proceeded with the examination of residual stress due to differential growth in adjoining tissue. It should be noted that the spherical symmetry is preserved in the analysis presented in the second and third chapters. The stability analysis was not performed for the spherical growth. This is beyond the scope of the presented 129 chapter. The computational model presented in the fifth chapter gave us an opportunity to examine the stress and strain fields emerging beyond the symmetry-breaking post-bifurcations regimes. The results show that the volumetric growth leads to stress and strain fields in a complex fashion, particularly in the sulcal floor and gyral crown. Thus, this work allows us to speculate on the possible broader connections of the presented results to physical and biological phenomena. Biological growth and remodeling change the mechanical state of the soft tissue and lead to residual fields. The experimental findings support the existence of residual stress fields in the soft tissues including the cortex. Due to the limited applicability of the cutting experiment to the cortex, the residual stress fields are not fully comprehended. This motivates extensive utilization of computational techniques to simulate the cortical growth and folding process. The main motivation of these numerical models is generally to understand the folding process in various aspects including cortical growth and folding mechanisms, developmental abnormalities, and the role of mechanical quantities on morphology. However, the mechanical state of the brain soft tissue following the growth and extensive folding period is still scarce. The presence of residual stress holds notable significance for several reasons. The residual stresses affect the elastic properties of soft tissue such as stiffness and cause heterogeneous and anisotropic internal stress fields. Thus, the residual stress fields can have a role in the microme- chanical state of the tissue that has a beneficiary or adverse effects. For example, the residual stress in arteries makes the stress distribution more uniform along the arterial wall [45]. Furthermore, the mechanical response of residually stressed tissues to applied loads (body forces and surface tractions) differs from that of stress-free tissues. Along with the effects on the mechanical state of the tissue, considering the role of mechanical properties of the tissues on the regulation of the essential behavior of the cellular structure, the residual fields have a potential role in mechanotransduction across all spatial scales of the brain. A local change in the residual stress field caused by an exogenous force below the threshold for mechanical failure may affect cellular physiology in the brain. The residual fields therefore should be included when seeking to model the micromechanical mechanisms that give rise to brain injury. 130 While the physical mechanisms of acute and secondary injuries led by extreme events (e.g., blunt impact, blast waves, cavitation, etc.) still remain unclear, it does seem clear that residual stresses could have a significant effect. Hence, the presence of residual fields in the reference state of the brain tissue should be regarded for the development of more elaborate models of the physical injury mechanisms. In that regard, further analysis is necessary for brain soft tissue possessing an intrinsic stress field associated with growth that is subject to a finite deformation (static or time-dependent). Such a study provides an understanding of the effect that residual stress has on the response of brain soft tissue to external loading that can lead to an injury. Similar to the framework presented in chapter 2, one possible approaches would be to develop a finite deformation boundary value problem for the geometry of mature tissue exhibiting residual stress field specifically caused by the morphogenesis. The reference state is actually deformed because of previous growth and folding processes, and intrinsically stressed even in the absence of body forces and surface tractions. It is therefore required to redefine reference configuration to account for the residual stress. To this end, one derives a general constitutive equation for a hyperelastic material with a residual stress field using the theory of invariants. The invariants are dependent upon the residual stress field, finite deformations by means of the right Cauchy- Green tensor, and the combination of them. There are 10 independent invariants in total for a compressible material in a general three-dimensional case, which reduces to 9 for the consideration of incompressible material. Cauchy stress is then expressed using a certain number of invariants involved to account for the residual stress adequately. An inherent challenge in this treatment lies in nonlinearity. It might not be possible to obtain an explicit analytical solution even if a simple form of strain energy function is utilized. The assumption of internal constraint alleviates the challenges to some extent at this point. However, considering the geometry of sulcus-gyrus formation and boundary conditions, the treatment might not be possible to solve analytically. Especially regarding our particular interest in intrasulcal deformations presented in Chapter 4, one can consider confining the attention to the intrasulcal 131 region instead of a relatively complex sulcus-gyrus formation. Later, the intrasulcal region can be idealized as illustrated in the figure 6.1. The intrasulcal region is here simplified as a cylindrical tube culminating in a hemispherical cap. In this representation, the cylindrical tube and hemispherical cap correspond to the sulcal walls and the depth of the sulcus, respectively. Figure 6.1 Idealized representation of the intrasulcal region. 6.1 On the elastodynamic analysis of incompressible hyperelastic bilayer spherical shell and its possible connections with residual fields In addition to the research work presented in this dissertation, we have conducted a study on the finite amplitude pure radial oscillation of multi-layered tubes. In this work, we have treated 132 White MatterInner RadiusOuter RadiusInterface RadiusGray MatterCerebrospinal FluidP(t) the elastodynamic problem of finite amplitude radial oscillation of a multilayer cylindrical tube where each layer consists of an incompressible hyperelastic multilayer. For the radial motion considered, the second-order equation of motion was derived. The combination of the symmetry of radial motion and the internal constraints of the material gives rise to a single degree of freedom system i.e., 𝑥 = 𝑥(𝑡) which is the oscillating location of the tube’s inner radius. The first integral of this equation yields the first-order differential equations - "energy equations", including the response integrals for each layer that embody the layer’s constitutive response through strain energy functions. This equation also enables the calculation of closed motion trajectories (orbits) in a phase plane with coordinates (𝑥, (cid:164)𝑥). Explicit results were then given for the free oscillations of a bilayer cylindrical tube with a Mooney-Rivlin type constitutive response for each layer to examine the influence of the layer differences for a specific problem. Specifically, we have investigated the effect of mismatch of shear modulus, mismatch of density, and interface location on the qualitative features of the phase plane orbits. One example of these numerical results is presented in Figure 6.2 which shows the motion as a trajectory in an (𝑥, 𝑣)-phase plane for a bilayer tube with a thick outer layer. In this specific case, the density and shear modulus of each layer are considered to be identical. To a large degree, common features were preserved through such parameter variation so long as the ratio of shear modulus to density is maintained among the various layers, which corresponds to the preservation of shear wave speed. The orbits were closed because the system is conservative. In the following part of the work, we have treated the elastodynamic analysis of the same motion using a direct energy method. The second-order radial equation of motion is thereby obtainable from either a Lagrangian or Hamiltonian dynamics formulation. The direct energy treatments offer a potentially more manageable alternative way to generate the specific elastic response integrals for each layer In addition, the energy formalism also provides a possible means to more easily incorporate additional phenomena such as residual stress. Further details about the treatments and numerical examples have been published in a research paper titled "Radially Oscillating Incompressible Hyperelastic Multi-Layer Tubes: Interface Effects and Energy Approach" [147]. 133 Figure 6.2 Demonstration of the motion as a trajectory in an (𝑥, 𝑣)-phase plane. In this example, the bilayer tube is considered with the dimensions for the inner radius 𝑅1 = 1, interface radius 𝑅2 = 2, and outer radius 𝑅3 = 20, The shear modulus 𝜇 and the density 𝜌 of each layer are considered to be identical i.e., 𝜇1 = 𝜇2 and 𝜌1 = 𝜌2. The table presents the turning point amplitudes (𝑎 and 𝑏), stroke (𝑏 − 𝑎) and the corresponding period of motion 𝜏. Energy formalism can be considered as an alternative treatment to exploring the possible effects of the residual fields on the response of brain soft tissue. For an illustrative explanation. Recalling the equation (3.23) in [147] for an illustrative explanation, the total stored energy 𝔘 for the problem at hand is expressed as 𝔘 = 𝔘ref + 𝜋𝑅4 1 𝑁 ∑︁ 𝑛=1 𝐻𝑛 (𝑥). (6.1.1) where 𝐻𝑛 (𝑥) is the base response function of layer 𝑛 with respect to non-dimensional inner radius stretch 𝑥, and 𝑅1 is the undeformed inner radius of the cylindrical tube, respectively. Here 𝔘ref denotes the baseline reference stored energy that can be utilized to quantify the residual stress in the material [63]. In case 𝔘ref = 0, the material does not possess any residual stress. It should be noted that the equation is derived for the problem described above, and presented to show how the effect of residual stress fields can be included in the analysis. Furthermore, It is important to remark that energy formalism allows one to obtain phase-plane trajectories for the considered motion. The 134 abstroke7.07680.04952.19732.1477τ aforementioned treatment would be more convenient for a stress analysis. The presented framework of radial oscillation of incompressible hyperelastic multilayer tubes can be easily implemented to examine specific examples associated with brain injury dynamics. Not only the free oscillations, but the problem can also be formulated for the forced oscillations with various pressure boundary conditions such as Heaviside step loading, blast loading, harmonic loading, and periodic step pulse load [114]. It is useful to remind that the symmetry of the motion provides simplification to the problem to some extent. Hence, the idealized intrasulcal region can be regarded for such an elastodynamic analysis, as given in figure 6.1. In this setting, the cylindrical tube represents the sulcal walls while the hemispherical end represents the depth of the sulcus. Assuming that the sulcal wall is a bilayer cylindrical tube with a thin layer (cortex) over the thick outer layer (subcortex) with similar material properties, figure 6.2 becomes a representative illustration of the elastodynamic response of a sulcal wall without any residual stress. 6.2 Potential roles of residual fields to brain injury mechanisms This section is reserved to speculate about the possible implementation of the mechanomor- phogenesis model presented in Chapter 5 to brain injury dynamics. Specifically, the potential roles of residual field induced by morphogenesis are discussed in two broad aspects: beneficial and adverse effects of intrinsic mechanical field on physical injury mechanisms seen in TBI, and mechanotransduction mechanisms. Potential roles on the physical injury mechanisms In fracture mechanics, the role of the maximum principal stress significantly influences crack behavior and its propensity to propagate. In the context of injury mechanics in soft tissues, the manifestation of cracks finds equivalence in the context of rips or tears resulting from the deformation. When the brain soft tissue is subject to a finite deformation, the internal mechanical state might have the potential to affect the occurrence of the injury independent of the loading condition itself. For instance, a preexisting tensile residual stress could accelerate the formation of micro-fissures during an episode of physical trauma, whereas a preexisting compressive residual stress field could provide some benefit in delaying fissure formation. The results indicated that the 135 maximum principal stresses in tension are localized around the gray-white matter interface and the fissure is likely to occur especially along the interfaces (see Figure 5.12). Parallel to the findings presented here, recent experimental studies suggest that extreme events such as cavitation [37, 22] and blast waves [121] lead to tissue rips and tears, particularly along the interfaces. However, there is insufficient data to substantiate the formation of tears and rips in relation to micro-fissures and the intrinsic mechanical state of the tissue, including residual stress and strains. Still, the extension of fracture mechanics to injury mechanics, incorporating residual fields, holds promise for elucidating the mechanisms of rips and tears observed in soft tissues subjected to extreme loads. Further exploration is warranted for the precise quantification of residual stress and strain fields, as well as the development of physical injury models that incorporate these intrinsic fields. Possible roles on mechanotransduction mechanisms Along with the role of residual fields on the continuum response of brain soft tissue, such a discussion on the role of intrinsic effects can be extended at the sub-tissue scale since the microenvironment plays a critical role in many cell-intrinsic mechanisms. The tissue comprises cells and an extracellular matrix (ECM), where the ECM serves as a scaffold offering mechanical support and driving biological signaling in cells within tissues. Recognized as a pivotal aspect of the microenvironment, the ECM regulates cell behaviors and phenotypes. Cells, in turn, possess the ability to sense and respond to the mechanical attributes of their microenvironment, known as mechanosensing. Through adaptive behaviors, cells govern critical aspects of cellular physiology, including membrane integrity, cellular morphology, and structural robustness. Experimental findings indicate that cells can recognize deformation energy resulting from matrix stretch, suggesting that cell mechanosensing is influenced by the deformation energy in the ECM, not solely its stiffness [94]. It is therefore reasonable to posit that the presence of growth-induced residual stress (strain) fields can be recognized by cell mechanosensing, and subsequently adapted into mechanochemical signaling pathways and cellular structures. This mechanical equilibrium between cells and residual fields thus constitutes a state of mechanical homeostasis, preserving structural integrity and functionality. However, the influence of mechanical 136 microenvironmental conditions in the presence of residual stress in the ECM remains unclear. Returning to the outcomes of the presented study in Chapter 5, the results indicate the existence of inhomogeneous mechanical intrinsic fields within gyrus-sulcus formations. This heterogeneity might suggest spatially varying homeostatic conditions for extra and intracellular structures along the gyrus-sulcus formation. Considering this diversity in cellular response, even acute physical forces below the threshold for mechanical failure could potentially induce local alterations in the microenvironment, disrupting mechanical homeostasis [60]. 137 BIBLIOGRAPHY [1] Erin L Abner, Peter T Nelson, Frederick A Schmitt, Steven R Browning, David W Fardo, Lijie Wan, Gregory A Jicha, Gregory E Cooper, Charles D Smith, Allison M Caban-Holt, et al. Self-reported head injury and risk of late-life impairment and ad pathology in an ad center cohort. Dementia and geriatric cognitive disorders, 37(5-6):294–306, 2014. [2] Milton Abramowitz and Irene A Stegun. Handbook of Mathematical Functions with Formu- las, Graphs, and Mathematical Tables. National Bureau of Standards Applied Mathematics Series 55. Tenth Printing. ERIC, 1972. [3] [4] A. Agosti, A. L. Gower, and P. Ciarletta. The constitutive relations of initially stressed incompressible mooney-rivlin materials. Mechanics Research Communications, 93:4–10, 2018. Howard G Allen. Analysis and design of structural sandwich panels. New York, NY: Pergamon Press, 1969. [5] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E Rognes, and Garth N Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015. [6] Martine Ben Amar and Alain Goriely. Growth and instability in elastic tissues. Journal of the Mechanics and Physics of Solids, 53(10):2284–2319, 2005. [7] [8] [9] Nickie N Andescavage, Adre Du Plessis, Robert McCarter, Ahmed Serag, Iordanis Evan- gelou, Gilbert Vezina, Richard Robertson, and Catherine Limperopoulos. Complex trajecto- ries of brain development in the healthy human fetus. Cerebral Cortex, 27(11):5274–5283, 2017. Donald H Barron. An experimental analysis of some factors involved in the development of the fissure pattern of the cerebral cortex. Journal of Experimental Zoology, 113(3):553–581, 1950. PV Bayly, RJ Okamoto, G Xu, Y Shi, and LA Taber. A cortical folding model incorporating stress-dependent growth explains gyral wavelengths and stress patterns in the developing brain. Physical biology, 10(1):016005, 2013. [10] PV Bayly, LA Taber, and CD Kroenke. Mechanical forces in cerebral cortical folding: a review of measurements and models. Journal of the mechanical behavior of biomedical materials, 29:568–581, 2014. [11] Nicola Bellomo, NK Li, and Ph K Maini. On the foundations of cancer modelling: se- lected topics, speculations, and perspectives. Mathematical Models and Methods in Applied Sciences, 18(04):593–646, 2008. [12] Maurice A Biot. Folding instability of a layered viscoelastic medium under compression. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 242(1231):444–454, 1957. 138 [13] Maurice A Biot. Surface instability of rubber in compression. Applied Scientific Research, Section A, 12(2):168–182, 1963. [14] Maurice A Biot. Edge buckling of a laminated medium. International Journal of Solids and Structures, 4(1):125–137, 1968. [15] Fabian Brau, Hugues Vandeparre, Abbas Sabbah, Christophe Poulard, Arezki Boudaoud, and Pascal Damman. Multiple-length-scale elastic instability mimics parametric resonance of nonlinear oscillators. Nature Physics, 7(1):56–60, 2011. [16] Nicholas J Braun, Katherine R Yao, Patrick W Alford, and Dezhi Liao. Mechanical injuries of neurons induce tau mislocalization to dendritic spines and tau-dependent synaptic dys- function. Proceedings of the National Academy of Sciences, 117(46):29069–29079, 2020. [17] Christopher Earls Brennen and Christopher E Brennen. Fundamentals of multiphase flow. Cambridge university press, 2005. [18] Silvia Budday, Ellen Kuhl, and John W Hutchinson. Period-doubling and period-tripling in growing bilayered systems. Philosophical Magazine, 95(28-30):3208–3224, 2015. [19] Silvia Budday, Charles Raybaud, and Ellen Kuhl. A mechanical model predicts morpholog- ical abnormalities in the developing human brain. Scientific reports, 4(1):1–7, 2014. [20] Silvia Budday and Paul Steinmann. On the influence of inhomogeneous stiffness and growth International Journal of Solids and on mechanical instabilities in the developing brain. Structures, 132:31–41, 2018. [21] Silvia Budday, Paul Steinmann, and Ellen Kuhl. The role of mechanics during brain devel- opment. Journal of the Mechanics and Physics of Solids, 72:75–92, 2014. [22] Saranya Canchi, Karen Kelly, Yu Hong, Michael A King, Ghatu Subhash, and Malisa Sarntinoranont. Controlled single bubble cavitation collapse results in jet-induced injury in brain tissue. Journal of the mechanical behavior of biomedical materials, 74:261–273, 2017. [23] Yanping Cao and John W. Hutchinson. Wrinkling Phenomena in Neo-Hookean Film/Substrate Bilayers. Journal of Applied Mechanics, 79(3):031019, 04 2012. [24] Yuli Cao, Mårten Risling, Elisabeth Malm, Anders Sondén, Magnus Frödin Bolling, and Mattias K Sköld. Cellular high-energy cavitation trauma–description of a novel in vitro trauma model in three different cell types. Frontiers in Neurology, 7:10, 2016. [25] Nancy Carney, Annette M Totten, Cindy O’Reilly, Jamie S Ullman, Gregory WJ Hawryluk, Michael J Bell, Susan L Bratton, Randall Chesnut, Odette A Harris, Niranjan Kissoon, et al. Guidelines for the management of severe traumatic brain injury. Neurosurgery, 80(1):6–15, 2017. [26] M. M. Carroll. Pressure maximum behavior inflation of incompressible elastic hollow spheres and cylinders. Quarterly of Applied Mathematics, 45:141–154, 1987. 139 [27] MM Carroll. Pressure maximum behavior in inflation of incompressible elastic hollow spheres and cylinders. Quarterly of applied mathematics, 45(1):141–154, 1987. [28] Miguel Cerrolaza, Sandra Shefelbine, and Diego Garzón-Alvarado. Numerical methods and advanced simulation in biomechanics and biological processes. Academic Press, 2017. [29] Xi Chen and John W Hutchinson. Herringbone buckling patterns of compressed thin films on compliant substrates. J. Appl. Mech., 71(5):597–603, 2004. [30] Yi-chao Chen and Anne Hoger. Constitutive functions of elastic materials in finite growth and deformation. In Advances in Continuum Mechanics and Thermodynamics of Material Behavior, pages 175–193. Springer, 2000. [31] Je G Chi, Elizabeth C Dooling, and Floyd H Gilles. Gyral development of the human brain. Annals of Neurology: Official Journal of the American Neurological Association and the Child Neurology Society, 1(1):86–93, 1977. [32] Cheng-Jen Chuong and Yuan-Cheng Fung. Residual stress in arteries. Frontiers in biome- chanics, pages 117–129, 1986. [33] Lucas da Costa Campos, Raphael Hornung, Gerhard Gompper, Jens Elgeti, and Svenja Caspers. The role of thickness inhomogeneities in hierarchical cortical folding. NeuroImage, 231:117779, 2021. [34] Mohsen Darayi, Mia E Hoffman, John Sayut, Shuolun Wang, Nagehan Demirci, Jack Con- solini, and Maria A Holland. Computational models of cortical folding: A review of common approaches. Journal of Biomechanics, page 110851, 2021. [35] Mohsen Darayi, Mia E Hoffman, John Sayut, Shuolun Wang, Nagehan Demirci, Jack Con- solini, and Maria A Holland. Computational models of cortical folding: A review of common approaches. Journal of Biomechanics, 139:110851, 2022. [36] Colette Dehay, Pascale Giroud, Michel Berland, Herbert Killackey, and Henry Kennedy. Contribution of thalamic input to the specification of cytoarchitectonic cortical fields in the primate: effects of bilateral enucleation in the fetal monkey on the boundaries, dimen- sions, and gyrification of striate and extrastriate cortex. Journal of Comparative Neurology, 367(1):70–89, 1996. [37] Carey E Dougan, Zhaoqiang Song, Hongbo Fu, Alfred J Crosby, Shengqiang Cai, and Shelly R Peyton. Cavitation induced fracture of intact brain tissue. Biophysical Journal, 121(14):2721–2729, 2022. [38] Laurence T Dunn. Raised intracranial pressure. Journal of Neurology, Neurosurgery & Psychiatry, 73(suppl 1):i23–i27, 2002. [39] Brian T Fagan, Sikhanda S Satapathy, J Neal Rutledge, and Steven E Kornguth. Simulation of the strain amplification in sulci due to blunt impact to the head. Frontiers in neurology, page 998, 2020. 140 [40] Benjamin Falcon, Jasenko Zivanov, Wenjuan Zhang, Alexey G Murzin, Holly J Garringer, Ruben Vidal, R Anthony Crowther, Kathy L Newell, Bernardino Ghetti, Michel Goedert, et al. Novel tau filament fold in chronic traumatic encephalopathy encloses hydrophobic molecules. Nature, 568(7752):420–423, 2019. [41] A. Font, N. K. Jha, H. Dehghani, J. Reinoso, and J. Meredio. Modeling of residually stressed, extended and inflated cylinders with applications to aneurysms. Mechanics Research Communications, 111:103643, 2021. [42] Giulia Franceschini, Davide Bigoni, Peter Regitnig, and Gerhard A Holzapfel. Brain tissue deforms similarly to filled elastomers and follows consolidation theory. Journal of the Mechanics and Physics of Solids, 54(12):2592–2620, 2006. [43] YC Fung and SQ Liu. Change of residual strains in arteries due to hypertrophy caused by aortic constriction. Circulation research, 65(5):1340–1349, 1989. [44] Yuan Cheng Fung. What are the residual stresses doing in our blood vessels? Annals of biomedical engineering, 19(3):237–249, 1991. [45] Yuan-cheng Fung. Biomechanics: mechanical properties of living tissues. Springer Science & Business Media, 2013. [46] Guangqiang Geng, Leigh A Johnston, Edwin Yan, Joanne M Britto, David W Smith, David W Walker, and Gary F Egan. Biomechanisms for modelling cerebral cortical folding. Medical image analysis, 13(6):920–930, 2009. [47] Mazdak Ghajari, Peter J Hellyer, and David J Sharp. Computational modelling of traumatic brain injury predicts the location of chronic traumatic encephalopathy pathology. Brain, 140(2):333–343, 2017. [48] Jacques Goeller, Andrew Wardlaw, Derrick Treichler, Joseph O’Bruba, and Greg Weiss. Investigation of cavitation as a possible damage mechanism in blast-induced traumatic brain injury. Journal of neurotrauma, 29(10):1970–1981, 2012. [49] A. Goriely and M. Ben Amar. On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity. Biomechanics and Modeling in Mechanobi- ology, 6:289–296, 2007. [50] Alain Goriely. The mathematics and mechanics of biological growth, volume 45. Springer, 2017. [51] Alain Goriely, Marc GD Geers, Gerhard A Holzapfel, Jayaratnam Jayamohan, Antoine Jérusalem, Sivabal Sivaloganathan, Waney Squier, Johannes AW van Dommelen, Sarah Waters, and Ellen Kuhl. Mechanics of the brain: perspectives, challenges, and opportunities. Biomechanics and modeling in mechanobiology, 14:931–965, 2015. [52] Kun Gou and Thomas J Pence. Hyperelastic modeling of swelling in fibrous soft tissue with application to tracheal angioedema. Journal of mathematical biology, 72(1):499–526, 2016. 141 [53] A. L. Gower, P. Ciarletta, and M. Destrade. Initial stress symmetry and its application in elasticity. Proc Roy Soc London, Ser A, 471(2183):20150448, 2015. [54] A. L. Gower, T. Shearer, and P. Ciarletta. A new restriction for initially stressed elastic solids. Quarterly Journal of Mechanics and Applied Mathematics, 40:455–478, 2017. [55] Albert Edward Green and RT Shield. Finite elastic deformation of incompressible isotropic bodies. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 202(1070):407–419, 1950. [56] S. E. Greenwald, Jr. Moore, J. E., A. Rachev, T. P. C. Kane, and J.-J. Meister. Experi- mental Investigation of the Distribution of Residual Strains in the Artery Wall. Journal of Biomechanical Engineering, 119(4):438–444, 11 1997. [57] Hans Gregersen, Ib Hansen. Morphometry and residual strain in porcine ureter. Scandinavian journal of urology and nephrology, 33(1):10–16, 1999. [58] Alexander Greiner, Stefan Kaessmair, and Silvia Budday. Physical aspects of cortical folding. Soft Matter, 17(5):1210–1222, 2021. [59] HC Han and YC Fung. Residual strains in porcine and canine trachea. Journal of biome- chanics, 24(5):307–315, 1991. [60] Matthew A Hemphill, Stephanie Dauth, Chung Jong Yu, Borna E Dabiri, and Kevin Kit Parker. Traumatic brain injury and the neuronal microenvironment: a potential role for neuropathological mechanotransduction. Neuron, 85(6):1177–1192, 2015. [61] J Claude Hemphill III, Steven M Greenberg, Craig S Anderson, Kyra Becker, Bernard R Bendok, Mary Cushman, Gordon L Fung, Joshua N Goldstein, R Loch Macdonald, Pamela H Mitchell, et al. Guidelines for the management of spontaneous intracerebral hemorrhage: a guideline for healthcare professionals from the american heart association/american stroke association. Stroke, 46(7):2032–2060, 2015. [62] Michel A Hofman. On the evolution and geometry of the brain in mammals. Progress in neurobiology, 32(2):137–158, 1989. [63] A. Hoger. On the residual stress possible in an elastic body with material symmetry. Archive for Rational Mechanics and Analysis, 88:271–290, 1985. [64] A. Hoger. On the determination of residual stress in an elastic body. J. Elasticity, 16:303–324, 1986. [65] Evan Hohlfeld and Lakshminarayanan Mahadevan. Unfolding the sulcus. Physical review letters, 106(10):105702, 2011. [66] Maria A Holland, Kyle E Miller, and Ellen Kuhl. Emerging brain morphologies from axonal elongation. Annals of biomedical engineering, 43:1640–1653, 2015. 142 [67] Ruoyu Huang, Raymond W Ogden, and Raimondo Penta. Mathematical modelling of residual-stress based volumetric growth in soft matter. Journal of Elasticity, 145(1):223– 241, 2021. [68] Ioan Humphreys, Rodger L Wood, Ceri J Phillips, and Steven Macey. The costs of traumatic brain injury: a literature review. ClinicoEconomics and outcomes research: CEOR, 5:281, 2013. [69] Khalid Iqbal, Cheng-Xin Gong, and Fei Liu. Hyperphosphorylation-induced tau oligomers. Frontiers in neurology, 4:112, 2013. [70] Fatemeh Jafarabadi, Shuolun Wang, and Maria A Holland. A numerical study on the influence of cerebrospinal fluid pressure on brain folding. Journal of Applied Mechanics, 90(7):071006, 2023. [71] Mir Jalil Razavi, Tuo Zhang, Tianming Liu, and Xianqiao Wang. Cortical folding pattern and its consistency induced by biological growth. Scientific reports, 5(1):14477, 2015. [72] Lihua Jin, Shengqiang Cai, and Zhigang Suo. Creases in soft tissues generated by growth. EPL (Europhysics Letters), 95(6):64002, 2011. [73] B. E. Johnson and A. Hoger. The dependence of the elasticity tensor on residual stress. Journal of Elasticity, 33:145–165, 1993. [74] Byron E Johnson and Anne Hoger. The dependence of the elasticity tensor on residual stress. Journal of Elasticity, 33:145–165, 1993. [75] Byron E Johnson and Anne Hoger. The use of a virtual configuration in formulating constitutive equations for residually stressed elastic materials. Journal of elasticity, 41:177– 215, 1995. [76] Joseph Kerwin, Atacan YÃŒcesoy, Suhas Vidhate, Bianca M. Dávila-Montero, Jacob L. Van Orman, Thomas J. Pence, Michaelann Tartis, Ricardo MejÃa-Alvarez, and Adam M. Willis. Sulcal cavitation in linear head acceleration: Possible correlation with chronic traumatic encephalopathy. Frontiers in Neurology, 13, 2022. [77] Vitaly A Klyachko and Charles F Stevens. Connectivity optimization and the positioning of cortical areas. Proceedings of the national academy of sciences, 100(13):7937–7941, 2003. [78] Steven Kornguth, Neal Rutledge, Gabe Perlaza, James Bray, and Allen Hardin. A proposed mechanism for development of cte following concussive events: head impact, water hammer injury, neurofilament release, and autoimmune processes. Brain sciences, 7(12):164, 2017. [79] Ji Lang, Rungun Nathan, Dong Zhou, Xuewei Zhang, Bo Li, and Qianhong Wu. Cavitation causes brain injury. Physics of Fluids, 33(3):031908, 2021. [80] PB Medawar Le Gros Clark. Deformation patterns on the cerebral cortex. In Essays on Growth and Form, pages 1–22, 1945. 143 [81] Bo Li, Yan-Ping Cao, Xi-Qiao Feng, and Huajian Gao. Mechanics of morphological insta- bilities and surface wrinkling in soft materials: a review. Soft Matter, 8(21):5728–5745, 2012. [82] Chih-Jen Lin and Jorge J Moré. Newton’s method for large bound-constrained optimization problems. SIAM Journal on Optimization, 9(4):1100–1127, 1999. [83] Eriko Sato Matsuo and Toyoichi Tanaka. Patterns in shrinking gels. Nature, 358(6386):482– 485, 1992. [84] Ann C McKee, Thor D Stein, Christopher J Nowinski, Robert A Stern, Daniel H Daneshvar, Victor E Alvarez, Hyo-Soon Lee, Garth Hall, Sydney M Wojtowicz, Christine M Baugh, et al. The spectrum of disease in chronic traumatic encephalopathy. Brain, 136(1):43–64, 2013. [85] [86] J. Merodio and R. W. Ogden. Extension, inflation and torsion of a residually-stressed circular cylindrical tube. Continuum Mechanics and Thermodynamics, 28:157–174, 2016. J. Merodio, R. W. Ogden, and J. Rodriguez. The influence of residual stress on finite deformation elastic response. International Journal of Non-Linear Mechanics, 56:43–49, 2013. [87] Julius Miklowitz. Elastic waves and waveguides. 1978. [88] Scott T Miller, Candice F Cooper, Paul Elsbernd, Joseph Kerwin, Ricardo Mejia-Alvarez, and Adam M Willis. Localizing clinical patterns of blast traumatic brain injury through computational modeling and simulation. Frontiers in Neurology, 12, 2021. [89] Philip H Montenigro, Charles Bernick, and Robert C Cantu. Clinical features of repetitive traumatic brain injury and chronic traumatic encephalopathy. Brain pathology, 25(3):304– 317, 2015. [90] Jingxin Nie, Lei Guo, Gang Li, Carlos Faraco, L Stephen Miller, and Tianming Liu. A computational model of cerebral cortex folding. Journal of theoretical biology, 264(2):467– 478, 2010. [91] Raymond W Ogden. Non-linear elastic deformations. Courier Corporation, 1997. [92] JH Omens, AD McCulloch, and JC Criscione. Complex distributions of residual stress and strain in the mouse left ventricle: experimental and theoretical models. Biomechanics and modeling in mechanobiology, 1(4):267–277, 2003. [93] Matthew B Panzer, Barry S Myers, Bruce P Capehart, and Cameron R Bass. Development of a finite element model for blast brain injury and the effects of csf cavitation. Annals of biomedical engineering, 40(7):1530–1544, 2012. [94] Valeria Panzetta, Sabato Fusco, and Paolo A Netti. Cell mechanosensing is regulated by substrate strain energy rather than stiffness. Proceedings of the National Academy of Sciences, 116(44):22004–22013, 2019. 144 [95] T. J. Pence and I. S. Wichman. Essential Mathematics for Engineers and Scientists. Cam- bridge University Press, 2020. [96] Thomas J Pence and Kun Gou. On compressible versions of the incompressible neo-hookean material. Mathematics and Mechanics of Solids, 20(2):157–182, 2015. [97] Thomas J Pence and Hungyu Tsai. Swelling-induced cavitation of elastic spheres. Mathe- matics and mechanics of solids, 11(5):527–551, 2006. [98] TJ PENCE and J SONG. Buckling instabilities in a thick elastic three-ply composite plate under thrust. International journal of solids and structures, 27(14):1809–1828, 1991. [99] WD Pilkey, DF Pilkey, and ZM Bi. Peterson?s stress concentration factors; 4th version, 2020. [100] Luka Pocivavsek, Robert Dellsy, Andrew Kern, Sebastián Johnson, Binhua Lin, Ka Yee C Lee, and Enrique Cerda. Stress and fold localization in thin elastic membranes. Science, 320(5878):912–916, 2008. [101] Pasko Rakic. Specification of cerebral cortical areas. Science, 241(4862):170–176, 1988. [102] Rea Ravin, Nicole Y Morgan, Paul S Blank, Nitay Ravin, Hugo Guerrero-Cazares, Alfredo Quinones-Hinojosa, and Joshua Zimmerberg. Response to blast-like shear stresses associated with mild blast-induced brain injury. Biophysical journal, 117(7):1167–1178, 2019. [103] Mir Jalil Razavi, Tuo Zhang, Xiao Li, Tianming Liu, and Xianqiao Wang. Role of mechanical factors in cortical folding development. Physical Review E, 92(3):032701, 2015. [104] Mir Jalil Razavi, Tuo Zhang, Tianming Liu, and Xianqiao Wang. Cortical folding pattern and its consistency induced by biological growth. Scientific reports, 5(1):1–14, 2015. [105] Isabel Reillo, Camino de Juan Romero, Miguel Ángel García-Cabezas, and Víctor Borrell. A role for intermediate radial glia in the tangential expansion of the mammalian cerebral cortex. Cerebral cortex, 21(7):1674–1694, 2011. [106] David P Richman, R Malcolm Stewart, John Hutchinson, and Verne S Caviness Jr. Mechan- ical model of brain convolutional development: Pathologic and experimental data suggest a model based on differential growth within the cerebral cortex. Science, 189(4196):18–21, 1975. [107] R. S. Rivlin and D. W. Saunders. Large elastic deformations of isotropic materials vii: Experiments on the deformation of rubber. Phil. Trans. R. Soc. Lond., A243:251–288, 1951. [108] E. K. Rodriguez, A. Hoger, and A. D. McCulloch. Stress-dependent finite growth in soft elastic tissues. Journal of Biomechanics, 27:455–467, 1994. [109] Edward K Rodriguez, Anne Hoger, and Andrew D McCulloch. Stress-dependent finite growth in soft elastic tissues. Journal of biomechanics, 27(4):455–467, 1994. 145 [110] J. Rodriguez and J. Merodio. Helical buckling and postbuckling of pre-stressed cylindrical tubes under finite torsion. Finite Elements in Analysis and Design, 112:1–10, 2016. [111] Robert S Salzar, Derrick Treichler, Andrew Wardlaw, Greg Weiss, and Jacques Goeller. Ex- perimental investigation of cavitation as a possible damage mechanism in blast-induced Journal of neurotrauma, traumatic brain injury in post-mortem human subject heads. 34(8):1589–1602, 2017. [112] Susan C Schwerin, Elizabeth B Hutchinson, Kryslaine L Radomski, Kapinga P Ngalula, Carlo M Pierpaoli, and Sharon L Juliano. Establishing the ferret as a gyrencephalic animal model of traumatic brain injury: optimization of controlled cortical impact procedures. Journal of neuroscience methods, 285:82–96, 2017. [113] Jeong-Sun Seo, Seungbok Lee, Jong-Yeon Shin, Yu Jin Hwang, Hyesun Cho, Seong-Keun Yoo, Yunha Kim, Sungsu Lim, Yun Kyung Kim, Eun Mi Hwang, et al. Transcriptome analy- ses of chronic traumatic encephalopathy show alterations in protein phosphatase expression associated with tauopathy. Experimental & molecular medicine, 49(5):e333–e333, 2017. [114] Moshen Shahinpoor and JL Nowinski. Exact solution to the problem of forced large am- plitude radial oscillations of a thin hyperelastic tube. International Journal of Non-Linear Mechanics, 6(2):193–207, 1971. [115] M. Shams, M. Destrade, and R. W. Ogden. Initial stresses in elastic solids: constitutive laws and acoustoelasticity. Wave motion, 48:552–567, 2011. [116] M. Shams and R. W. Ogden. On rayleigh type surface waves in an initially stressed incom- pressible elastic solid. IMA Journal of Applied Mathematics, 79:360–376, 2014. [117] Richard Skalak. Growth as a finite displacement field. In Proceedings of the IUTAM symposium on finite elasticity, pages 347–355. Springer, 1981. [118] Richard Skalak, G Dasgupta, M Moss, E Otten, P Dullemeijer, and H Vilmann. Analytical description of growth. Journal of theoretical biology, 94(3):555–577, 1982. [119] IH Smart and GM McSherry. Gyrus formation in the cerebral cortex of the ferret. ii. description of the internal histological changes. Journal of anatomy, 147:27, 1986. [120] G Elliot Smith. The central nervous system in: Cunningham’s textbook of anatomy. by Arthur Robinson. New York, 1928. [121] Miguel A Gama Sosa, Rita De Gasperi, Alejandro J Paulino, Paul E Pricop, Michael C Shaughness, Eric Maudlin-Jeronimo, Aaron A Hall, William GM Janssen, Frank J Yuk, Nathan P Dorr, et al. Blast overpressure induces shear-related injuries in the brain of rats exposed to a mild traumatic brain injury. Acta neuropathologica communications, 1:1–15, 2013. [122] Georg F Striedter, Shyam Srinivasan, and Edwin S Monuki. Cortical folding: when, where, how, and why? Annual review of neuroscience, 38:291–307, 2015. 146 [123] Ghatu Subhash, Saranya Canchi, Yu Hong, Malisa Sarntinoranont, and Michael A King. Damage in brain tissue due to single bubble cavitation shock. In Mechanics of Biological Systems and Materials, Volume 6, pages 1–5. Springer, 2016. [124] Dhananjay Radhakrishnan Subramaniam, Ginu Unnikrishnan, Aravind Sundaramurthy, Jose Enrique Rubio, Vivek Bhaskar Kote, and Jaques Reifman. Cerebral vasculature in- fluences blast-induced biomechanical responses of human brain tissue. Frontiers in bioengi- neering and biotechnology, page 1057, 2021. [125] Tuomas Tallinen, Jun Young Chung, John S Biggins, and L Mahadevan. Gyrification from constrained cortical expansion. Proceedings of the National Academy of Sciences, 111(35):12667–12672, 2014. [126] Tuomas Tallinen, Jun Young Chung, François Rousseau, Nadine Girard, Julien Lefèvre, and Lakshminarayanan Mahadevan. On the growth and form of cortical convolutions. Nature Physics, 12(6):588–593, 2016. [127] XG Tan, AJ Przekwas, and RK Gupta. Computational modeling of blast wave interaction with a human body and assessment of traumatic brain injury. Shock Waves, 27(6):889–904, 2017. [128] PH Todd. A geometric model for the cortical folding pattern of simple folded brains. Journal of theoretical biology, 97(3):529–538, 1982. [129] Roberto Toro and Yves Burnod. A morphogenetic model for the development of cortical convolutions. Cerebral cortex, 15(12):1900–1913, 2005. [130] Timothy J Van Dyke and Anne Hoger. A new method for predicting the opening angle for soft tissues. J. Biomech. Eng., 124(4):347–354, 2002. [131] David C Van Essen. A tension-based theory of morphogenesis and compact wiring in the central nervous system. Nature, 385(6614):313–318, 1997. [132] Theodore B VanItallie. Traumatic brain injury (tbi) in collision sports: Possible mechanisms of transformation into chronic traumatic encephalopathy (cte). Metabolism, 100:153943, 2019. [133] SN Verner and K Garikipati. A computational study of the mechanisms of growth-driven folding patterns on shells, with application to the developing brain. Extreme Mechanics Letters, 18:58–69, 2018. [134] J Vossoughi. Ii. intimal residual stress and strain in large arteries. In 1993 Bioengng. Conf. ASME, 1993. [135] Linlin Wang, Jianyao Yao, and Ning Hu. A mechanical method of cerebral cortical folding development based on thermal expansion. Scientific Reports, 9(1):1914, 2019. [136] Qiming Wang and Xuanhe Zhao. A three-dimensional phase diagram of growth-induced surface instabilities. Bulletin of the American Physical Society, 60, 2015. 147 [137] Qiming Wang and Xuanhe Zhao. Beyond wrinkles: Multimodal surface instabilities for multifunctional patterning. Mrs Bulletin, 41(2):115–122, 2016. [138] Shuolun Wang, Nagehan Demirci, and Maria A Holland. Numerical investigation of biome- chanically coupled growth in cortical folding. Biomechanics and Modeling in Mechanobi- ology, pages 1–13, 2020. [139] Shuolun Wang, Nagehan Demirci, and Maria A Holland. Numerical investigation of biome- chanically coupled growth in cortical folding. Biomechanics and Modeling in Mechanobi- ology, 20:555–567, 2021. [140] Xiaoyu Wang, Julien Lefèvre, Amine Bohi, Mariam Al Harrach, Mickael Dinomais, and François Rousseau. The influence of biophysical parameters in a biomechanical model of cortical folding patterns. Scientific Reports, 11(1):1–14, 2021. [141] Johannes Weickenmeier, Rijk de Rooij, Silvia Budday, Paul Steinmann, Timothy C Ovaert, and Ellen Kuhl. Brain stiffness increases with myelin content. Acta biomaterialia, 42:265– 272, 2016. [142] Anna Wermer, Joseph Kerwin, Kelsea Welsh, Ricardo Mejia-Alvarez, Michaelann Tartis, and Adam Willis. Materials characterization of cranial simulants for blast-induced traumatic brain injury. Military medicine, 185(Supplement_1):205–213, 2020. [143] Gang Xu, Philip V Bayly, and Larry A Taber. Residual stress in the adult mouse brain. Biomechanics and modeling in mechanobiology, 8:253–262, 2009. [144] Gang Xu, Andrew K Knutsen, Krikor Dikranian, Christopher D Kroenke, Philip V Bayly, and Larry A Taber. Axons pull on the brain, but tension does not drive cortical folding. Journal of biomechanical engineering, 132(7), 2010. [145] Atacan Yucesoy. The role of morphology and residual stress on blast-induced traumatic brain injury. Master’s thesis, Michigan State University. Mechanical Engineering, 2022. [146] Atacan Yucesoy and Thomas J Pence. On the inflation of residually stressed spherical shells. Journal of Elasticity, pages 1–20, 2021. [147] Atacan Yucesoy and Thomas J Pence. Radially oscillating incompressible hyperelastic multi-layer tubes: Interface effects and energy approach. Journal of Elasticity, pages 1–20, 2023. [148] Qiu Yue, Kim Sangwoo, and Thomas J Pence. Plane strain buckling and wrinkling of neo-hookean laminates. International journal of solids and structures, 31(8):1149–1178, 1994. [149] V. Zamani and T. J. Pence. Swelling, inflation, and a swelling-burst instability in hyperelastic spherical shells. International Journal of Solids and Structures, 125:134–149, 2017. [150] Vahid Zamani and Thomas J Pence. Swelling, inflation, and a swelling-burst instability in hyperelastic spherical shells. International Journal of Solids and Structures, 125:134–149, 2017. 148 [151] Jan Zavodnik, Andrej Košmrlj, and Miha Brojan. Rate-dependent evolution of wrinkling films due to growth on semi-infinite planar viscoelastic substrates. Journal of the Mechanics and Physics of Solids, 173:105219, 2023. [152] Tuo Zhang, Mir Jalil Razavi, Xiao Li, Hanbo Chen, Tianming Liu, and Xianqiao Wang. Mechanism of consistent gyrus formation: an experimental and computational study. Scien- tific reports, 6(1):1–11, 2016. [153] Feng Zhu, Christina Wagner, Alessandra Dal Cengio Leonardi, Xin Jin, Pamela VandeVord, Clifford Chou, King H Yang, and Albert I King. Using a gel/plastic surrogate to study the biomechanical response of the head under air shock loading: a combined experimental and numerical investigation. Biomechanics and modeling in mechanobiology, 11(3):341–353, 2012. [154] KA Zimmerman, J Kim, C Karton, L Lochhead, DJ Sharp, T Hoshizaki, and M Ghajari. Player position in american football influences the magnitude of mechanical strains produced in the location of chronic traumatic encephalopathy pathology: A computational modelling study. Journal of biomechanics, 118:110256, 2021. 149