Influence of Initial Residual Stress on Growth and Pattern Creation for a Layered Aorta

Residual stress is ubiquitous and indispensable in most biological and artificial materials, where it sustains and optimizes many biological and functional mechanisms. The theory of volume growth, starting from a stress-free initial state, is widely used to explain the creation and evolution of growth-induced residual stress and the resulting changes in shape, and to model how growing bio-tissues such as arteries and solid tumors develop a strategy of pattern creation according to geometrical and material parameters. This modelling provides promising avenues for designing and directing some appropriate morphology of a given tissue or organ and achieve some targeted biomedical function. In this paper, we rely on a modified, augmented theory to reveal how we can obtain growth-induced residual stress and pattern evolution of a layered artery by starting from an existing, non-zero initial residual stress state. We use experimentally determined residual stress distributions of aged bi-layered human aortas and quantify their influence by a magnitude factor. Our results show that initial residual stress has a more significant impact on residual stress accumulation and the subsequent evolution of patterns than geometry and material parameters. Additionally, we provide an essential explanation for growth-induced patterns driven by differential growth coupled to an initial residual stress. Finally, we show that initial residual stress is a readily available way to control growth-induced pattern creation for tissues and thus may provide a promising inspiration for biomedical engineering.

www.nature.com/scientificreports www.nature.com/scientificreports/ for bio-tissues. As shown by Fig. 1(A,B), the residual stress in biological tissues can be released only partially by cutting them in different directions. A complete release of residual stress would require an infinite number of cuts.
Similarly, some investigations on stress-dependent or strain-dependent growth law cannot be set up properly without using the actual residual stress. By dividing the total growth deformation into many small steps that are further decomposed by the MD model, Goriely and Ben Amar 22 proposed an incremental, cumulative, and continuous growth law, allowing some growth steps to develop from a residually stressed state. With the exception of the very first step, every step grows from a residually stressed state, which provides a way to analyze the influence of residual stress on the growth process. However, the residual stress in any accumulative growth step stems from the previous growth step; most importantly, it is still necessary to assume an initial residual stress-free state in the first step, something which, again, is very difficult to prescribe and measure in real biological tissues.
It follows that the MD growth model must be modified to account for the existing initial residual stresses found in living, growing continuum tissues, which result from complicated growth and remodelling processes or other bio-interactions 6,23 .
Hence, taking the initial state to be an arbitrary residual stress state and building on the residual stress theory of Hoger et al. 21,[24][25][26] , we recently proposed 27 a modified multiplicative decomposition growth (MMDG) model. In this paper we set out to show the influence of initial residual stress on growth-induced residual stress and the following pattern creation by focusing on a scenario modeling a three-dimensional, layered, initially stressed, growing aorta.
We adopt some existing experimental and theoretical results of aortic residual stress distributions for our initial residual stress, and discuss its influence by way of a magnitude factor. We compute the distributions of growth-induced residual stress and the geometrical changes of the aorta starting from different initial residual stress states. Based on a subsequent stability analysis for the growing aorta with initial residual stress, we obtain the critical wrinkling solution using linear incremental theory and surface impedance methods. We investigate the influence of the initial residual stress level and of the differential growth extent on pattern development for the layered aorta and provide sound explanations for their effect on pattern creation.
Finally, we emphasise that our results may provide an effective strategy to obtain targeted patterns by controlling the initial residual stress instead of the differential growth or swelling extent 19,20 . Hence by going beyond bio-tissue growth, this work can also provide inspiration for industrial manufacturing, nano-fabrication, or self-assembly of soft solids which display behaviors similar to growth or swelling [28][29][30][31][32] .

MMDG Framework and Basic Equations
Based on the concept of multiplicative decomposition of the total growth process, the overall growth process for an initially stressed material is also decomposed into separate parts. As shown in Fig. 2, our innovation is to introduce an initial elastic deformation F 0 for releasing the initial residual stress to a virtual stress-free configuration. Thereafter, the pure growth deformation F g can happen between two virtual stress-free and incompatible configurations. Finally, a pure elastic deformation F e is induced when the tissue restores the compatibility of the bio-tissue after growth, which is the direct source of the associated residual stress.
According to this new description of the growth process, the total growth deformation is decomposed as Here we assume that the material in the virtual stress-free state is a natural material following the constitutive equation of an incompressible neo-Hookean solid, also adhering to the initial stress symmetry condition 26  4) were almost stress-free whilst the core layers (2 and 3) were subjected to large inhomogeneous axial and radial residual stresses; (C) Layer-specific structure of a cut human aorta ring (from Holzapfel et al. 1 ), revealed once it is unloaded from blood pressure and from axial and radial residual stresses: the circumferential residual stress causes a consequent contraction in the unloaded state, leading to buckling and even delimitation on the inner surface.
www.nature.com/scientificreports www.nature.com/scientificreports/ although other models could equally have been chosen, e.g. 25,[33][34][35] . Hence, from Eq. (1) follows the constitutive equation for the Cauchy stress σ of growing bio-tissues with initial residual stress τ as where p, p 0 are the Lagrange multipliers in the current and reference configurations, respectively. Here we assume isotropic growth, so that = = where N and n are the normal vectors to the surfaces in the reference and current configurations, respectively.

A Growing Aorta with An Experimentally-Determined Initial Residual Stress Field
Layered tubular tissues are the focus of this paper, as they are common in the body, e.g. artery, airway, intestine, etc. They are all known to possess a high degree of residual stress, as shown experimentally by cuts. It is also known that their bio-functions depend significantly on their residual stress levels and their morphology. Figure 1(C) shows a ring of an aged human aorta, which can be modelled as a bi-layered tube according to Holzapfel and collaborators 1,2 .
Here we follow their model and take the aorta as a bi-layer, by combining the intima and the media into one inner layer and having the adventitia as the outer layer in a three-dimensional Euclidean space. The reference configuration with initial residual stress is associated with the cylindrical coordinates Θ R Z ( , , ), and the current configuration with a new residual stress state is associated with the coordinates (r, θ, z).
Holzapfel and collaborators 1,2 determined the distribution of residual stress in aged human aortas from measurements of the circumferential opening angles and axial bending angles for each layer, see Fig. 3(A) for the protocol. We reproduced their results (details not shown here) to get τ in non-dimensional form, see Fig. 3(B) (and compare with the dimensional version 2 ). From here on, we multiply this distribution by a factor α to introduce a magnitude factor of the initial residual stress and clearly illustrate its effects: α = 0 corresponds to the stress-free state, α = 1 corresponds to the initial stress of Holzapfel and collaborators 1,2 . We take their shear moduli μ = .
34 4 kPa in for the inner layer and μ = . 17 3 kPa ad for the adventitia layer, respectively. In addition, for convenience and better display of the results, we scale down the dimensions of the bi-layered aorta by a factor 10, resulting in the following geometry: = .
R 0 5911 mm for the inner, interface, and outer radii, respectively.
The radial stress is continuous across the aortic thickness and remains compressive throughout, while the circumferential and axial stresses are discontinuous at the interface between the layers. The circumferential stress changes from compressive on the inner surface of the intima to tensile at the interface, and remains tensile in the adventitial layer. The axial stress is compressive in the inner layer and tensile in the adventitial layer. The axial stress has the largest magnitude, and the circumferential stress is also consequent; the radial stress is almost negligible in comparison, see Fig. 3

(B).
According to the above bi-layered modeling, we now impose the constrained growth deformation (overall growth deformation) gradient tensors in the residually stressed inner and adventitia layers by prescribing www.nature.com/scientificreports www.nature.com/scientificreports/ , where k = in, ad denotes the corresponding quantities in the inner layer and adventitia, respectively. Then, based on Eq. (2), we compute the following non-zero Cauchy stress components in the current state, where p k 0 , p k are yet undetermined variables: p k depends on the boundary conditions in the current configuration, and p k 0 is determined from the incompressibility condition for the natural materials, Now that we have the distribution of the initial residual stress, p k 0 is found in terms of τ k and shear moduli μ k as 3, 3 k , and τ I 1, k, τ I 2, k, τ I 3, k are the principal invariants of the initial residual stress tensor. Then, from the incompressibility condition we find that   in adventitia layer, and different initial residual stress magnitudes α = 1 (same as Holzapfel and collaborators 1,2 ) and α = 2 (twice their level). We see that the magnitude of the residual stress components increases with increasing growth volume in the inner layer. Moreover, the changes of residual stresses due to growth are very similar between the two different initial residual stress levels α = 1 and α = 2, which demonstrates that the level of self-balanced initial residual stress does not affect the accumulation law of growth-induced residual stress. However, owing to the superposition effect of initial residual stress and growth-induced residual stress, there will be a significant difference in the total residual stress.  but has some slight impact on the growth process by the little difference of the absolute radii. This difference could also play a role in creating wrinkles and patterns, see the next section.
We expect that increasing compressible circumferential stress and axial stress on the inner surface eventually induces a buckling of the surface and creates wrinkles or folds, a phenomenon that is commonly found in most tubular tissues. However, we see from Fig. 4 that the circumferential tensile stress in the inner layer (intima and media) and the tensile axial stress in the outer layer (adventitia) are increased, which could have either opposite effects on instability or help the tissues keep stable.
Based on the MD growth model, Ciarletta et al. 20 have shown that geometrical dimensions and stiffness contrast between the layers play a significant role in growth-induced pattern selection for tubular tissues. The results above have demonstrated the influence of the initial residual stress on growth-induced residual stress accumulation and on geometry changes. It follows that the initial residual stress must have a certain impact on the subsequent pattern creation, and we now turn to wrinkling analysis to determine its exact influence.

Directional Pattern Creation by initial Residual Stress and Differential Growth
Morphogenesis and pattern formation are closely related to several specific bio-functions. It is thus imperative to study the relationships between growth, initial stress, change in shape and eventually, pattern selection. Some pioneering works studied pattern selection for tubular tissues, where the tubular organ grows from a stress-free initial state. Hence Ciarletta et al. 20 showed that increasing the thickness or stiffness ratio between the outer and inner tubular layers creates fewer wrinkles and folds in the circumferential direction but more wrinkles in the axial direction, thus providing an insight into bio-medical engineering applications to achieve directional control or selection of growth-induced pattern creation. Here by setting α = 0, we expect to recover their results. Our innovation is a proposal to directionally control or select patterns by initial residual stress distribution (α ≠ 0), instead of altering the geometry or elasticity of the tissues.
We describe the wrinkles by a three-dimensional incremental field θ θ θ = + + where m and κ π = n 2 / are the circumferential and longitudinal wave numbers, respectively ( is the current length of the tube). In other words, m and n are the numbers of circumferential and axial wrinkles, respectively. Then the amplitude functions U, V, W, Q are determined by solving numerically the incremental equilibrium equations and boundary conditions. This task is best achieved by using the Stroh formulation, which we recall in the appendix. Then, by iterating over the wave numbers m and k, we eventually get the critical value of each case which relates to the growth-induced pattern creation.
First we need to evaluate the critical level of initial residual stress signalling when the aorta is unstable in its reference configuration before any differential growth takes place (for other examples of this situation, see 23,36,37 ). Here the initial residual stress is based on experimental data and we quantify its influence by the magnitude factor , in the presence of an initial residual stress found in an aged human aorta 1,2 , with magnitudes α = 1 and α = 2; (B) Changes in radii and thicknesses of the bi-layered aorta with increasing inner layer growth and different initial residual stress magnitudes α = 1 and α = 2.
www.nature.com/scientificreports www.nature.com/scientificreports/ α. Since there are no wrinkling patterns of note in an aged human aorta, we must check that the initial residual state of Holzapfel et al. 1,2 at α = .
1 0 is in the stability range for α. Figure 5(A) displays how the critical magnitude value α cr relates to the possible wrinkling modes, with number of circumferential wrinkles = … m 0, 1, , 4 and number of axial wrinkles = … n 1, 2, , 6. The figure shows that the smallest α cr for any given number of axial wrinkles n occurs when there are no circumferential wrinkles ( = m 0): hence based only on an increase of initial residual stress (and no growth), the bi-layered aorta would buckle in the axial direction. Moreover, we also find the minimal critical value α = . min( ) 2 343 cr , happening with = m 0, = n 5. It means that the surfaces of the initially residually stressed aorta (at α = . 1 0) are smooth and free of wrinkles according to this specific modeling, in line with experimental observation.
From now on we limit the span of magnitude α for the initial residual stress to be less than min(α cr ), to avoid wrinkles in the initial configuration. Figure 5(B) then shows the influence of α on growth-induced pattern creation. As α increases toward min(α cr ), the amount of differential growth (volume accumulation in the inner layer) decreases, tending to 1 (no growth) eventually, as expected. In that limit we recover = m 0 cr , = n 5 cr of Fig. 5(A). The theoretical and numerical analysis in the appendix reveals that critical pattern modes strongly depend on both the magnitude α of the initial residual stress and the subsequent differential growth, rather than only on the latter as in previous studies. Figure 5 shows that growth-induced wrinkle creation is significantly altered by changing slightly the magnitude of initial residual stress and then letting the aorta grow at a constant rate. With no growth ( = J 0 g in ), the initially stressed aorta can only wrinkle axially. With no initial residual stress (α = 0), differential growth can only lead to a = = m n 1 pattern. By combining both effects, we can see circumferential and axial wrinkles arise and add constructively to form a further 2D pattern with = m 2, = n 1, see Fig. 6(A) for examples. Figure 6(B) shows a comparison of the cross-section areas of the inner volume between an atherosclerotic and a healthy aorta 38 , which is qualitatively consistent with the conclusion drawn from Fig. 6(A) when the initial residual stress gets bigger. So our results may provide another possible explanation of physical pathology for atherosclerosis.

Concluding Remarks
We demonstrate the influence of initial residual stress on the growth and pattern creation of layered aortas. A modified multiplicative decomposition growth (MMDG) framework is employed to incorporate and consider initial residual stress when we simulate the growth process theoretically. We highlight the differences in residual stress accumulation and in geometrical changes due to the presence of initial residual stress. We also implement a linearized stability analysis to unveil pattern creation at critical levels of growth and initial stress.
As a conclusion, we point out that understanding the effect of the initial residual stress provides an alternative way to select objective patterns or to control some bio-functions via morphological change. This may also present an inspiring insight for directional bionic self-assembly or robot manufacturing by initial residual stress.

Method
We did not conduct ourselves the experiments on human aortas reported here. They were conducted by the authors of reference 1 in 2007, and we are simply reporting some of their findings and using them as a base to derive some of our results. The use of autopsy material from human subjects in reference 1 was approved in 2007 by the Ethics Committee, Medical University Graz, Austria, see details in reference 1 .

Appendix: Wrinkling Analysis
Here, we use an infinitesimal elastic deformation χ′ as a perturbation of the large deformation solution in the current configuration. The corresponding incremental displacement gradients  F (with respect to the reference configuration) and  F I (with respect to the current configuration) are related as: =   F F F I . By assuming that the incremental deformation is infinitesimal and transient, and that the corresponding growth process is independent of the stress and strain fields, we have the relationship: =   F F F e I e for the pure growth gradient.
The incremental nominal stress component in push-forward form is are the instantaneous elastic moduli 25 . The non-zero incremental equilibrium equations are along the three principal directions:    To solve these equations numerically, subject to the appropriate incremental boundary conditions, we use the surface impedance method, see details elsewhere [39][40][41] .