Homogenized constrained mixture models for anisotropic volumetric growth and remodeling

Constrained mixture models for soft tissue growth and remodeling have attracted increasing attention over the last decade. They can capture the effects of the simultaneous presence of multiple constituents that are continuously deposited and degraded at in general different rates, which is important to understand essential features of living soft tissues that cannot be captured by simple kinematic growth models. Recently the novel concept of homogenized constrained mixture models was introduced. It was shown that these models produce results which are very similar (and in certain limit cases even identical) to the ones of constrained mixture models based on multi-network theory. At the same time, the computational cost and complexity of homogenized constrained mixture models are much lower. This paper discusses the theory and implementation of homogenized constrained mixture models for anisotropic volumetric growth and remodeling in three dimensions. Previous constrained mixture models of volumetric growth in three dimensions were limited to the special case of isotropic growth. By numerical examples, comparison with experimental data and a theoretical discussion, we demonstrate that there is some evidence raising doubts whether isotropic growth models are appropriate to represent growth and remodeling of soft tissue in the vasculature. Anisotropic constrained mixture models, as introduced in this paper for the first time, may be required to avoid unphysiological results in simulations of vascular growth and remodeling.


Introduction
Growth and remodeling of biological soft tissues have attracted increasing interest over the last two decades. For a comprehensive recent review, the reader is referred to (Cyron and Humphrey 2016). Naturally associated with large deformations, they are best addressed within the mathematical framework of nonlinear continuum mechanics. To this end, Rodriguez et al. developed a kinematic growth theory based on a simple multiplicative split of the deformation gradient into an elastic part and an inelastic (growth) part, similar to mathematical models of plasticity or viscoelastic fluids (Rodriguez et al. 1994). Thanks to its conceptual simplicity and computational efficiency, this kinematic growth theory is still widely used (Ambrosi et al. 2011;Menzel and Kuhl 2012) and has been applied to various types of soft tissues and biomaterials (e.g., Göktepe 2010;Zöllner 2012Zöllner , 2013Ambrosi and Guana 2007;Grytz 2012;Albero et al. 2014;Ben Amar and Goriely 2005;Goriely and Vandiver 2010;Böl and Bolea Albero 2014;Taber and Eggers 1996;Ambrosi and Mollica 2004;Vandiver and Goriely 2009). It suffers, however, from a couple of severe limitations. In particular, it neglects the simultaneous production and degradation processes of multiple constituents (such as elastin and collagen) that is characteristic of living tissues. These processes are key to understand and capture phenomena such as the change of the opening angle of arteries after application of proteases (Cyron et al. 2016a) or in hypertension (Matsumoto and Hayashi 1996). To incorporate these processes, so-called constrained mixture models were developed (Humphrey and Rajagopal 2002), motivated by the idea of multi-network theory (Rajagopal and Wineman 1992). These models have now been established as a second major approach to growth and remodeling in soft tissues (Figueroa 2009;Baek et al. 2006;Valentin 2009;Wilson et al. 2013;Zeinali-Davarani et al. 2011;Holzapfel 2008, 2007). They describe growth and remodeling as a superposition of continuous degradation and deposition processes of differential mass increments of multiple constituents in each differential volume element. To avoid the high complexity and computational cost of these models, Watton et al. developed, motivated by Humphrey (1999), an alternative constrained mixture model that is not based on multi-network theory but only tracks for each constituent a so-called recruitment stretch that accounts for deposition and degradation of mass increments in a lumped manner (e.g., Eriksson 2014; Watton et al. 2004Watton et al. , 2011. With simple kinematic growth models, both isotropic and anisotropic volumetric growth can be modeled straightforwardly. For constrained mixture models, the implementation of volumetric growth is more challenging. Such models have therefore most often been combined with mechanical membrane models, where growth was simply described as a change of the thickness parameter of the membrane (Figueroa 2009;Baek et al. 2006;Wilson et al. 2013;Watton et al. 2004Watton et al. , 2011Wilson et al. 2012;Zeinali-Davarani and Baek 2012;Watton and Hill 2009). After some extensions to thick-walled axisymmetric blood vessels (Virag 2015;Karšaj et al. 2010), only recently constrained mixture models for volumetric growth in general geometries were proposed (Eriksson 2014;Grytsan et al. 2015;Valentín et al. 2013). They were, however, limited to isotropic volumetric growth.
Recently the so-called homogenized constrained mixture models were introduced (Cyron et al. 2016b). They are essentially based on the same micromechanical model of growth and remodeling as classical constrained mixture models based on multi-network theory. However, using a form of temporal homogenization across differential mass increments deposited at different times, they can capture the effects of growth and remodeling by a split into an elastic and inelastic deformation gradient of each constituent comparable to the one used in kinematic growth models. This way, homogenized constrained mixture models combine the realistic micromechanical foundation of classical constrained mixture models based on multi-network theory with the conceptual simplicity and low computational cost of kinematic growth models. Moreover, they enable a straightforward incorporation of both isotropic and anisotropic volumetric growth. By Cyron et al. (2016b) mainly the mathematical and mechanical foundations of homogenized constrained mixture models were introduced, illustrated only by the simple example of growth in a two-dimensional curved membrane continuum.
In this paper, we discuss in detail how homogenized constrained mixture models can be applied to study anisotropic volumetric growth in three dimensions. We demonstrate general good agreement but also some quantitative differences between the results of three-dimensional solid models compared to equivalent two-dimensional membrane models. This supports on the one hand the appropriateness of previous membrane models for general academic studies but also indicates that three-dimensional models may be an indispensable tool to ensure the model accuracy that may be required in future clinical applications. Moreover, we present some evidence raising doubts whether isotropic volumetric growth models are appropriate to study vascular growth and remodeling in three dimensions. Isotropic models apparently exhibit some difficulties in reproducing the growth response of arteries to hypertension, and they also appear to suffer from an inherent susceptibility to mechanobiological instability.
In this paper, first-and second-order tensors are printed in boldface and double contraction products of two tensors are denoted by a colon and time derivatives by an overdot. For constrained mixtures composed of n constituents (i.e., material types like collagen or elastin), a superscript i = 1, . . . , n denotes a mechanical quantity of the ith constituent. Determinants of second-order tensors are denoted by |.|, traces by tr (.). For tensor products, the symbol ⊗ is used.

Kinematics and balance equations
Let B 0 ⊂ R 3 be a reference configuration of a continuous body mapped at time t ≥ 0 to a current configuration B(t). Each material point X ∈ B 0 is mapped to its current position x (t, X) ∈ B(t) by with deformation gradient Without loss of generality, we assume herein the initial condition B(0) = B 0 . Note that B 0 can be an arbitrary configuration and need not satisfy any specific mechanical conditions. In practice, B 0 will be most often the in vivo configuration of a biological tissue before some specific event starts a certain growth and remodeling process of interest. For example, it can be the healthy configuration of a blood vessel right before the vessel suffers some mechanical or biochemical damage that initiates the growth of an aneurysm. Reference volume elements dV in B 0 are mapped to volume elements dv = J dV in the current configuration B(t) with the Jacobian determinant J = |F| > 0. Each volume element consists of a constrained mixture of n different constituents (e.g., collagen, elastin, or smooth muscle) that are supposed to undergo the same deformation over time. That is, the total deformation gradient F in (2) is the same for all constituents. The local stress-free state may, however, differ between the constituents and in fact in general even between the differential mass increments of these constituents deposited at different times τ . That is, for each differential mass increment of the ith constituent deposited at time τ we split the deformation gradient into an individual elastic part F i(τ ) e and inelastic part F Here gr captures the differences between the local stressfree configurations of different mass increments due to deposition at a different time and in a different configuration as well as continuous growth and remodeling. Growth and remodeling happen on slow time scales (at least hours, typically days to months) and can thus be modeled as a quasi-static process so that dynamic effects like inertia or viscoelasticity can usually be neglected and the balance of linear momentum becomes in the interior of the body div (P) with body force vector b 0 (per unit mass) and the first Piola-Kirchhoff stress P. The total mass density (per unit reference volume) is the sum of the mass densities i 0 of the individual constituents. Following Eqs. (2.51) and (4.6) from Holzapfel (2000), the mass density in the current spatial configuration is computed from the mass density in the reference configuration by Note that Equation (4.6) in Holzapfel (2000) was derived under the assumption of conservation of mass but that it holds true also generally if the amount of mass in reference volume elements changes over time so that it can be used herein. The first Piola-Kirchhoff stress P = ∂Ψ/∂F can be computed from the total strain energy function (per unit reference volume) where Ψ i and W i are the strain energies per unit reference volume and mass of each constituent, respectively. W i depends in general only on the elastic deformation of each mass increment, that is, on the elastic right Cauchy-Green tensor or, equivalently, on its invariants. For isotropic materials, there exist only three invariants, and for anisotropic materials, additional pseudo-invariants have to be considered (cf. section 6 in Holzapfel 2000). Below, we will use not only the first Piola-Kirchhoff stress P but also the second Piola-Kirchhoff stress where E = F T F − I /2 is the Green-Lagrange strain tensor, I is the identity tensor, and is the second Piola-Kirchhoff stress of the ith constituent and ϕ i its volume fraction in the current configuration. The push-forward of S i renders which can be interpreted as the Cauchy stress in the fictitious case that the volume fraction ϕ i of the ith constituent equals one. In a given configuration B(t), the balance of linear momentum (4) can be evaluated at any point in time using gr are known or if at least some additional equations are given to compute them at any time. The above equations represent a common basis of all constrained mixture models for growth and remodeling proposed so far. The difference between these models is mainly the way how the i 0 and F i(τ ) gr are computed, respectively. In the following section, we will point out in detail how these quantities can be computed in homogenized constrained mixture models for anisotropic volumetric growth and remodeling.

General
In constrained mixture models, growth and remodeling are conceptualized as a simultaneous continuous degradation and deposition of n different constituents, which is also referred to as mass turnover. New mass is assumed to be deposited into the extant matrix with a certain prestress σ i pre at a rate˙ i 0+ ≥ 0. At the same time, extant mass is continuously degraded at a rate˙ i 0− ≥ 0 so that the net rate of mass production iṡ Given that in general the prestress σ i pre with which new mass is deposited is unequal to the current stress σ i of the extant mass that is being degraded, mass turnover changes in general the current average stress and thus also traction-free state of a constituent even when mass deposition and degradation balance (i.e.,˙ i 0 = 0). This change of the traction-free state in the absence of any change of local mass density or volume is necessarily associated with some change of the microstructure of the tissue that is referred to herein as remodeling. The local traction-free configuration of a constituent will in general not only change by this kind of turnover-based remodeling but also by growth, that is, a change of mass. If˙ i 0+ and˙ i 0− do not balance in (12) so that˙ i 0 = 0 the resulting change of mass is in general associated with a local change of the volume required to accommodate the mass in a certain region of the body. This growth-related local change of volume can be, similarly as the above-described remodeling, expressed as inelastic change of the local traction-free configuration of a certain constituent.
To account for changes of the traction-free state resulting from growth and remodeling, classical constrained mixture models keep track of all mass increments deposited at different times τ . For each constituent, the strain energy Ψ i is computed on the basis of a history integral incorporating each of these mass increments and its individual stress-free configuration. Computing this history integral is involved and time consuming. Therefore, Cyron et al. (2016b) introduced a homogenization across the different mass increments deposited at different times. This homogenization is based on the idea that the gross effect of growth and remodeling can be captured by some average inelastic deformation gradient F i gr for all mass increments of the ith constituent. That is, F i(τ ) gr = F i gr for all τ in (3) and we can rewrite (3) as where F i gr can be decomposed multiplicatively as Here the remodeling part F i r captures changes in the microstructure by mass turnover. The growth part F i g is directly related to a change of the mass per unit reference volume i 0 . Then, to close the system of equations (1)-(8) and enable a computation of the mechanical equilibrium configuration of the body at each time, one needs, in addition to proper initial and boundary conditions, two kinds of constitutive equations. First, one needs mechanical constitutive equations defining the strain energies W i of the constituents. Numerous such constitutive equations have been developed over the last two decades for soft tissues. All these equations could be used in the framework developed herein. For simplicity, we will rely in this paper largely on the ideas introduced by Holzapfel et al. (2000) and discuss their application briefly in Sect. 3.3 below. Second, one needs mechanobiological constitutive equations defining the evolution of the inelastic deformation gradients F i r and F i g as well as of the mass densities i 0 depending on mechanical stimuli such as stress or strain. These mechanobiological evolution equations will be discussed in the following Sect. 3.2.

Mass production
The most important mechanobiological constitutive equations are the ones defining the net mass production rates of each constituent in the current mechanical state. Following Cyron et al. (2016b) (and thereby the assumptions previously made, for example, also by Figueroa 2009;Wilson et al. 2013Wilson et al. , 2012, we assume that the net mass production (i.e., the difference between mass production and mass degradation) is governed by the difference between the co-rotated Cauchy stress tensor σ i R = R T σ i R and some (typically constant) homeostatic stress σ h . Here R is the orthonormal rotation tensor in the polar decomposition of the deformation gradient F = RU with U the material stretch tensor. Then the net mass production can be written aṡ with some (typically constant) gain-type second-order tensor K i σ . Here the first summand on the right-hand side is the mechano-regulated net mass production rate (depending on the current stress). The second summandḊ i (t) is a generic rate function that can be used to describe additional deposition or also damage processes that affect the net mass production but are not driven by Cauchy stress but other factors (such as mechanical fatigue and chemical degradation processes). For the theory developed herein also any other constitutive equation defining˙ i 0 could be chosen. Collagen and smooth muscle, the major constituents subject to growth and remodeling in load-bearing soft tissues such as arteries, are often modeled using quasi-one-dimensional fiber families. Then, in both the co-rotated Cauchy stress tensor σ i R and the homeostatic stress tensor σ i h only the components in reference fiber direction a i 0 are unequal to zero and denoted by σ i = a i 0 ⊗ a i 0 : σ i R and σ i h = a i 0 ⊗ a i 0 : σ i h , respectively, and one can rewrite (15) aṡ

Remodeling
For a known net mass production rate˙ i 0 and a deposition prestress σ i pre , it was demonstrated in Cyron et al. (2016b) that the evolution of the inelastic remodeling deformation gradient F i r is given by (cf. Eq. (17) and Eq. (31) in Cyron et al. (2016b) and note the symmetry of Here L i r =Ḟ i r (F i r ) −1 is the remodeling velocity gradient, S i = |F| F −1 σ i F −T and S i pre = |F| F −1 σ i pre F −T are the second Piola-Kirchhoff stress and deposition prestress, and T i is the average turnover time, that is, the period within which a mass increment is degraded and replaced by a new mass increment. This time is related to the half-life time measured experimentally, for example, by radioactive labeling by Nissen et al. (1978) and Etminan (2014), by the factor ln (2). Modeling mass turnover as a continuous replacement of fibers by new fibers in the same material direction, the rotational part of F i r is zero and F i r thus symmetric so that (17) provides in d dimensions a system of d (d + 1) /2 independent equations for computing the same number of unknown components ofḞ i r in any given mechanical state. S i pre (t) is the elastic stress of new mass that is produced at time t at a rate˙ i 0+ . At the same time, extant mass with stress S i is degraded at a rate˙ i 0− . Note that in our model we do not consider single mass increments produced and degraded but only keep track of the net change of the elastic stress in each constituent that results from such a combined production and degradation process. This net change of elastic stress can be expressed in terms of an inelastic remodeling rateḞ i r that is computed by (17). As shown in Proposition 1 of , the prestress σ i pre has to correspond to the homeostatic stress introduced in Sect. 3.2.1, that is,

Growth
An increase or decrease in mass of a differential material element in the body is in general associated with an increase or decrease of its volume, noting that the space a material element occupies depends (for a given density) on its mass. This change of volume can be captured by an inelastic growth deformation gradient F i g describing the change of size and shape of the differential volume element by deposition or degradation of mass. The local inelastic deformation gradient resulting from a certain change of mass depends in general on the micromechanical and geometric characteristics of the underlying growth process. The development of appropriate constitutive equations defining the evolution of F i g in time is thus non-trivial and a comprehensive discussion of the relation between the growth process on the micro-scale and the resulting F i g on the continuum-scale would go beyond the scope of this paper. Instead, we focus on a particularly simple model. We assume that deposition or degradation of mass of any constituent can be imagined as a process that changes (per unit time) the shape of volume elements in which it happens by an inelastic deformation gradient rateĠ i . This change of shape is assumed herein to occur by way of an inelastic reorganization of the whole volume element that affects all constituents equally. Therefore, all constituents are experiencing the same inelastic growth deformation gradient F g = F i g . Its rate equals the sum of the growth-related deformation gradient ratesĠ i contributed by the individual constituents: For living soft tissue, it is often assumed that growth occurs as a change of volume while the spatial density remains constant in time. Moreover, due to its high volume fraction of water, living soft tissue is typically modeled as an incompressible material, that is, its elastic deformation is supposed to be isochoric. This paper mainly focuses on the discussion of a theoretical framework rather than physiological details. Therefore, although more complex models could be thought of and might describe more realistically the relation between fluid and solid constituents in soft tissue, we assume herein for simplicity that also the elastic deformation of each individual constituent is isochoric, that is, F i e = 1. Isochoric elastic deformation implies via (17) automatically also an isochoric remodeling deformation gradient, that is, F i r = 1. Then, with (13) and (14) From (6), (19), and |F(0)| = 1, we conclude and thuṡ This equation has to hold for any net mass production rate of any constituent and thus also for the special case that only the net mass production of a single constituent whose index is without loss of generality i is not equal to zero. Then with (5),˙ 0 (t) =˙ i 0 (t) and all summands on the right-hand side of (21) except for the ith one are equal to zero so that (with Eq. (1.241) in Holzapfel 2000) i The rateĠ i can in general be written aṡ where the second-order tensor B i defines the growth direction and is normalized without loss of generality such that tr B i = 1. From (22) and (23), we conclude and thus, with (23) and (20) Note that (18) and (25) uniquely determine the evolution of the growth deformation gradient F i g for given mass density 0 (t), mass production rates˙ i 0 (t), and growth directions B i . The mass production rates (and thus also current mass) are determined by (15). Thus, only the growth direction B i has to be defined in addition to (15) or (16). Like (17) with regard to remodeling, (18) and (25) account for the effect of growth in a homogenized sense. That is, they do not consider separately the individual effects of all the mass increments deposited at different times in the tissue but only their gross effect by way of the total increase of mass.

Mechanical constitutive equations
Numerous mechanical constitutive equations for soft biological tissue have been published over the last decades. The homogenized constrained mixture model developed herein can be combined with arbitrary strain energy functions for the different constituents of the constrained mixture. Nevertheless, we will specify in this section a set of such functions that has been turned out particularly suitable for vascular biomechanics in the past (Wilson et al. 2013;Zeinali-Davarani et al. 2011;Wilson et al. 2012;) and which is therefore also used in the examples in Sect. 4. The arterial wall is modeled as a constrained mixture of elastin, (circumferential) smooth muscle and four families of collagen fibers oriented in circumferential, axial, and diagonal (+/ − 45 • ) direction (cf. Fig. 1). If in the following a certain parameter or quantity relates specifically to elastin, collagen, or smooth muscle, we will use the superscript el, co, or sm, respectively. In this section, we omit superscripts in F gr and F e , given that in each equation it is evident to which constituent in the constrained mixture they refer.
Collagen and the passive elasticity of smooth muscle are represented by quasi-one-dimensional fiber families, respectively. Assume these are aligned with the unit direction vector a 0 in reference configuration and the unit direction vector a gr = F gr a 0 / F gr a 0 in the inelastically deformed intermediate configuration, where F gr is the inelastic deformation gradient of the respective fiber family. Then we introduce the pseudo-invariant where C e is the elastic Cauchy-Green deformation tensor and λ e the elastic stretch in fiber direction of the respective constituent. The strain energy of a collagen fiber family is then a Fung-type exponential function with material parameters k co 1 and k co 2 . For smooth muscle, we use the strain energy function with a passive Fung-type part with material parameters k sm 1 and k sm 2 . Given the general importance of active smooth muscle tone (Murtada et al. 2012;Murtada and Holzapfel 2014;Murtada et al. 2010aMurtada et al. , b, 2015 in the vasculature, we also incorporate a simple model of active smooth muscle tension based on the strain energy function with λ act the active stretch in fiber direction, σ actmax the maximal active Cauchy stress, and λ m and λ 0 the active stretches at maximum and zero active stress. Consistent with other studies (Wilson et al. 2013), one assumes that λ act = λ/λ act whereλ act equals, due to fast muscle remodeling, at every point in time the total stretch in muscle direction in mechanical equilibrium so that in any equilibrium configuration with total stretch λ in fiber direction (compared to the reference configuration) we have λ act = 1 but yet ∂λ act /∂λ = 1/λ act = 1/λ, which is important to compute the stress resulting from (30). As pointed out already in Sect. 3.2.3, we assume herein that each individual constituent is incompressible under elastic deformation, that is, |F e | = 1. For constituents represented by uniaxial fiber families, (26)-(30) define via the equilibrium of linear momentum only the elastic stretch λ e in fiber direction. Incompressibility is then enforced by the simple kinematic assumption that the elastic stretch of the fiber material perpendicular to the fiber axis equals 1/ √ λ e so that the elastic deformation gradient of a fiber material becomes with identity tensor I.
In two-dimensional membrane models of arteries (e.g., Wilson et al. 2013;Zeinali-Davarani et al. 2011;Wilson et al. 2012;, elastin is often represented by an isotropic two-dimensional neo-Hookean strain energy function where C e is the elastic Cauchy-Green deformation tensor of elastin, a ⊥ gr defines a unit vector in wall thickness direction in the inelastically deformed intermediate configuration of elastin, A ⊥ gr = a ⊥ gr ⊗ a ⊥ gr is the structural tensor in wall thickness direction, and A gr = I − a ⊥ gr ⊗ a ⊥ gr a structural tensor for the tangential plane to the vessel wall. Note that a ⊥ gr = F gr a ⊥ 0 / F gr a ⊥ 0 with the unit vector a ⊥ 0 in wall thickness direction in the reference configuration and the inelastic deformation gradient F gr (capturing growth and remodeling) of elastin. Note that (32) is identical to, for example, Equation (2.6) in Wilson et al. (2013), that is, it exactly represents the strain energy that is typically used in two-dimensional membrane models for elastin, using only a notation that is generall enough to deal with arbitrary three-dimensional geometries. In this paper, we do not primarily focus on constitutive modeling of the arterial wall. Rather we aim at establishing a new mathematical framework for growth and remodeling in three dimensions. To test this framework, it is instructive to compare the results of our three-dimensional simulations with previously reported results from two-dimensional membrane models. Therefore, we herein do not fully exploit the freedom our three-dimensional model is granting to develop a detailed constitutive model that accounts for the inhomogeneity of the arterial wall in thickness direction. Instead we rather use a constitutive model as similar as possible to the ones used in previous two-dimensional models.
The above constitutive equations (26)-(32) are identical to the constitutive framework used previously in two dimensions (e.g., Wilson et al. 2013Wilson et al. , 2012Cyron et al. 2016b. The only modification we add to this framework in this paper is an additional strain energy contribution for elastin that ensures that the degrees of freedom that appear additionally in three dimensions compared to two dimensions are endowed with a proper stiffness (in order to avoid numerical problems). To this end, a three-dimensional isotropic neo-Hookean strain energy W el 3D is added to the elastin strain energy from (32) to constrain shear deformation. Thereby, one uses the isochoric elastic deformation gradient with F e the elastic deformation gradient of elastin and F e = 1. The related isochoric elastic right Cauchy-Green tensor of elastin isC e =F T eF e , and its first invariant is tr C e . Using these quantities, we define a three-dimensional isotropic neo-Hookean strain energy contribution of elastin in the form with the material parameter μ el 3D governing the shear stiffness. For elastin, incompressibility cannot be enforced by a simple kinematic assumption like (31). Therefore, we use the volumetric penalty function where a sufficiently high penalty parameter ε is used to enforce a nearly incompressible behavior of elastin. The total strain energy of elastin is then composed by the contributions from (32), (34), and (35) and becomes For a numerical implementation, the first and second derivatives of the above strain energy functions with respect to the Green-Lagrange strain E are required. These are provided in Appendix 2.

General
We implemented the homogenized constrained mixture model developed in Sect. 3 in our in-house research code BACI (written in C++) and studied volumetric growth and remodeling in a simple model aorta that is represented by a thick-walled cylinder of length L, inner radius R, and wall thickness H with an internal (blood) pressure p. This thick-walled cylinder was discretized with 8-noded hexahedral finite elements. An Fbar method was applied to avoid problems with locking (de Souza Neto et al. 1996). Based on a convergence study that was performed to ensure a negligible numerical error, we discretized the geometry with 552,960 elements (480 elements in axial direction, 6 elements in radial direction, and 192 elements in circumferential direction). The evolution of growth and remodeling (i.e., (16)-(18), and (25)) was solved using a backward Euler time integration scheme (cf. Appendix 3) with a time step size ten times smaller than the smallest turnover time T i (cf. (17)) of any constituent. The direction vectors of fiber families representing collagen and smooth muscle were averaged within each finite element. So far, volumetric constrained mixture models of arterial growth assumed isotropic growth (Eriksson 2014;Grytsan et al. 2015;Valentín et al. 2013).
However, experiments suggest in particular for arteries rather anisotropic growth, that is, material deposition or degradation (at least mainly) in thickness direction (Matsumoto and Hayashi 1996). Therefore, we study below both the case of isotropic growth represented by B i = I/3 in (25) and the case of anisotropic growth in thickness direction represented by B i = a ⊥ 0 ⊗ a ⊥ 0 with the unit vector a ⊥ 0 in wall thickness direction (in reference configuration). In case of isotropic growth, we conclude from (18) and (25) that the time derivative of the growth deformation gradient and thus, because of the initial condition F i g (t = 0) = I, also the growth deformation gradient itself can be written as some scalar prefactor times the identity tensor. With (20) this prefactor can be computed as 0 (t) / 0 (0) such that Similarly, one concludes in case of anisotropic growth in wall thickness direction There is some evidence that in adult arterial tissue mechanically functional elastin is not produced any longer but only slowly degraded with a half-life time of several decades (Cyron and Humphrey 2016). Thus, elastin growth can be characterized by the mass production rate˙ el 0 (t) = − el 0 (t) /T el rendering in (17) the remodeling velocity gradient L el r = 0 and thus F el r (t) = F el r (0) at any time. Collagen and smooth muscle are modeled by quasi-onedimensional fiber families. Note that in this case (17) defines the remodeling deformation gradient uniquely only via the additional kinematic assumption of incompressibility, that is, (31), which confines the space of possible elastic and remodeling deformations (cf. also Appendix 1 in Cyron et al. (2016b)).
Herein we simulate growth and remodeling starting from a homeostatic initial configuration (i.e., a mechanical equilibrium configuration where no mechano-regulated growth and remodeling occur). This means that for constituents subject to growth and remodeling (i.e., collagen and smooth muscle), the initial elastic stress or equivalently stretch has to equal the respective homeostatic one (which itself typically equals the deposition stretch, cf. Proposition 1 in ). It is then an important question how the elastic stretch of elastin has to be chosen so that under this constraint mechanical equilibrium in a given initial configuration is fulfilled. In a membrane, the answer to this question is straightforward. However, in a general three-dimensional solid it is in general non-trivial (cf. also page 10 in Eriksson (2014)). To overcome this problem, we use herein the procedure described in Appendix 1. Essentially, in this procedure we slightly vary the two-dimensional stiffness parameter μ el 2D and initial radial elastic stretch of elastin λ el⊥ e (t = 0) over the wall thickness such that both mechanical equilibrium and a homeostatic state of collagen and smooth muscle are ensured in the given initial configuration. Note that only these two parameters are varied in wall thickness direction, whereas all other parameters are kept constant to keep the model particularly similar to previously published membrane models and allow thereby instructive comparisons. All simulation parameters are summarized in Table 1 and chosen identical to the parameter values used in previous computational studies of aortic aneurysms with two-dimensional membrane models and reported in Table 1 of Wilson et al. (2013). Originally, these parameter values were determined by Wilson et al. (2012) from biaxial mechanical testing data reported by Geest et al. (2004) by nonlinear regression. The only additional parameters we used in this paper are μ el 2D and λ el⊥ e (t = 0), which are computed as pointed out in Appendix 1, and the damage parameters. The latter have been chosen heuristically for our specific computational example in Sect. 4.2 in order to resemble a typical damage in the vascular wall as it may occur during the initiation of an aneurysm.

Enlargement of axisymmetric model aneurysm
In the first example, we study the enlargement of an idealized fusiform abdominal aortic aneurysm over 15 years, starting from the above-described axisymmetric cylindrical geometry. Generally, anisotropic growth in thickness direction is assumed in this section [cf. (38) in Sect. 4.1]. Elastin is assumed to be subject to a damage rate functionḊ where we use a coordinate system with origin in the center of the cylinder and X 3 -axis aligned with the symmetry axis of the cylinder so that X 3 [−L/2; L/2]. The first summand on the right-hand side of (39) corresponds to a normal age-related degradation (driven by a Poisson process with time constant T el ), the second summand to an additional damage (e.g., by proteolytic activity in the vessel wall) that starts at t = 0 and spreads out along the vessel axis (beginning in the center) following a Gaussian-type function. It finally degrades up to a fraction D max of the initial elastin in the vessel wall. We compared the response of our threedimensional model aorta to the damage from (39) with the response of a two-dimensional membrane model where also a homogenized constrained mixture model for growth and remodeling was used (cf. section 3.2 of Cyron et al. (2016b)). We performed the comparison for a range of gain parameters k co σ = k sm σ ∈ [0.05, 0.07, 0.09, 0.11, 0.13, 0.15] /T co . The results are shown in Fig. 2. Generally, the agreement between the solution of the two-dimensional membrane model and the three-dimensional model developed herein is good. To examine the origin of the minor differences between both, we performed further computational studies. These suggested that the differences mainly result from the different ways how elastin is modeled in two and three dimensions, respectively. In the two-dimensional membrane model, the strain energy of elastin is given by the two-dimensional neo-Hookean function in (32) only where the elasticity parameter is chosen μ el 2D = 72 J/kg. In the three-dimensional model, the strain energy of elastin includes additionally a three-dimensional neo-Hookean contribution and a volumetric penalty term as can be seen in (36). Moreover, in three dimensions, we use μ el 3D = 72 J/kg and compute μ el 2D from the condition of mechanical equilibrium in the initial configuration, which provides values in the range 9 J/kg ≤ μ el 2D ≤ 40 J/kg. These minor differences between the model in two and three dimensions could not be avoided due to the additional degrees of freedom in three dimensions and the requirement to ensure initial mechanical equilibrium with respect to all of them. Yet, altogether the behavior of the two-and threedimensional model is very similar. In both models, the vessel responds to the elastin damage specified in (39) for small gain parameters by an unbounded enlargement (mechanobiological instability). By contrast, for large gain parameters, the blood vessel can apparently, after a transient period of enlargement, reattain a new stable and only slightly dilated configuration, which is in good agreement with the theory of mechanobiological stability proposed by Cyron et al. To distinguish between different collagen fiber families, we use, when appropriate, superscripts co(0 • ), co(90 • ), co(45 • ), and co(−45 • ) for circumferential, axial and diagonal fibers. For collagen and smooth muscle, the elastic stretch in the initial configuration equals the homeostatic stretch in fiber direction (that is equivalent to a certain homeostatic Cauchy stress). For μ el 2D and λ el⊥ e (t = 0) not a single value is given but rather a parameter interval within which these parameters are varied in wall thickness direction according to the procedure described in Appendix 1 in order to ensure both mechanical and mechanobiological equilibrium in the given initial geometry (2014), . To illustrate the simulations discussed in this section, the von Mises stress and the normalized collagen mass density at different points in time and for two different growth parameters k co σ = k sm σ = 0.15/T co and k co σ = k sm σ = 0.05/T co are shown in Figs. 3 and 4.
In previous three-dimensional constrained mixture models of arterial growth and remodeling (Eriksson 2014; Valen-  Grytsan et al. 2015), a loss of elastin was in general associated with a contraction of the vessel radius rather than a dilatation as one would intuitively expect after a weakening of the wall. The reason for this unphysical behavior was the way how volumetric growth was conceptualized in these models. The growth direction was implicitly determined by the incompressibility constraint using a volumetric deviatoric split of the deformation gradient. This way only isotropic growth could be modeled. So, if elastin was degraded and the total tissue volume decreased, the tissue contracted equally in all spatial directions and thus also in circumferential direction so that the radius signifi-cantly decreased. To overcome this apparently unphysical behavior, Eriksson (2014) introduced the distinction between constant individual density (CID) and adaptive individual density (AID) growth. In AID growth, a loss of mass does not change the tissue volume itself and so a loss of elastin does not entail a contraction of the artery. While this idea can, from a practical perspective, help avoid an unphysical contraction of the artery after a loss of elastin, its micromechanical foundation appears controversial in several respects. The medial layer of arteries is organized in a sandwich-like structure as discussed in detail by O'Connell (2008). Thin lamellar sheets of elastin spanning in circumferential and Void (that is, fluid-filled) spaces (white) are expected to rapidly fill by a relative motion of the adjacent material layers due to the pressure in radial direction throughout the wall, which leads to a reduced wall thickness as illustrated in c axial direction separate thick layers of smooth muscle and collagen fibers (cf. Fig. 5). A sudden loss of elastin mass would mean that the lamellar sheets are partially replaced by void (that is, only fluid-filled) spaces or at least spaces of massively reduced resistance against compression. Noting the pressure in the thickness direction of the vessel wall, these spaces would be expected to close rapidly by a relative motion of the adjacent material layers in wall thickness direction. Thus, it appears at least questionable whether a loss of elastin would really change the mass density of the tissue as implicitly assumed by the AID concept proposed in Eriksson (2014) or whether it would not rather locally shrink the tissue in thickness direction of the vessel wall. Such shrinking in thickness direction is exactly what the anisotropic volumetric growth model developed herein describes. It is important to note that previous isotropic volumetric constrained mixture growth models were in principle not able to capture a thinning of the vessel wall by a partial loss of tissue mass without altering significantly the traction-free configuration of the remaining mass in circumferential direction, which leads to the above-mentioned unphysical contraction. By contrast, the anisotropic growth model developed herein overcomes this problem in a fairly natural way (cf. Fig. 5). Matsumoto and Hayashi (1996) examined the response of the aorta to a sudden and permanent increase in blood pressure.

Thickening of arterial wall in hypertension
In this section, we study a similar situation with the parameters from Table 1. Unlike in Sect. 4.2, no damage of elastin is considered noting the short period studied compared to the long half-life time of elastin and the absence of particular Fig. 6 Adaptation of a healthy aorta to a minor increase in blood pressure δp = 20 mmHg for different gain parameters k co σ governing collagen production. The maximal inner radius is plotted over time. Note that the system is initially in a homeostatic state so that the maximal radius does not change until the system is perturbed at time t = 0. a Anisotropic growth: simulated vessels (solid lines) stabilize for k co σ = 0.34/T co around radius expected from experiments (dashed line). b While for anisotropic growth (dashed lines) vessel geometry stabilizes quickly after perturbation, isotropic growth (solid lines) leads for the same gain parameters k co σ = 0.30/T co , 0.65/T co , 1.00/T co to unstable growth and remodeling (continued dilatation of the vessel over time instead of stabilization around a new equilibrium state); this observation suggests a generally higher susceptibility of isotropic growth to mechanobiological instability proteolytic insults that could be expected to accelerate elastin degradation. At time t = 0, the blood vessel is subjected to a sudden increase in blood pressure p from 100mmHg to 120 mmHg. With a gain parameter k co σ = k sm σ = 0.34/T co and assuming anisotropic growth in thickness direction of the vessel wall, the vessel is observed to stabilize, after a transient growth period, in a new slightly dilated state that corresponds well to the one expected from the experimental observations of Matsumoto and Hayashi (1996) (cf. Fig. 6a). For a discussion how the data from experiments with rats in Matsumoto and Hayashi (1996) are compared with the dilatation of a human aorta as described by the parameters in Table 1, the reader is referred to section 3.3 of . Again the results obtained with the three-dimensional model developed in this paper are nearly identical to the ones obtained with the membrane model used by .
Interestingly, the results change completely if in the threedimensional model isotropic volumetric growth is assumed (similar to Eriksson 2014;Valentín et al. 2013;Grytsan et al. 2015). In this case, the vessel undergoes an unstable dilatation even for gain parameters as large as k co σ = k sm σ = 1.00/T co (cf. Fig. 6b). Not only is this behavior completely different from the stabilizing behavior observed experimentally by Matsumoto and Hayashi (1996), it also appears highly unphysiological. Blood vessels are expected to be subject to minor changes of blood pressure in numerous situations and if these always resulted in an unbounded unstable dilatation, aneurysms would be expected in nearly the whole and also younger population (rather than only 5-10% of the elderly population). The reason for the unstable dilatation in case of isotropic growth can be understood easily from the theory of mechanobiological stability . In this theory, one considers in general some state of mechanobiological equilibrium like the initial configuration in this example (where mechanical equilibrium is satisfied and no mechano-regulated growth and remodeling occur). Such an equilibrium state is called mechanobiologically stable if the system considered can reattain after small perturbations a mechanobiological equilibrium state identical to or close to the original one. In cylindrical blood vessels, after an increase in blood pressure, thickening of the wall helps restore a homeostatic stress state because it decreases the wall stress back to the initial level. On the other hand, circumferential dilatation of the wall further increases the wall stress due to Laplace' formula (which says that the wall stress increases linearly with the inner radius of the vessel). In case of anisotropic growth in thickness direction, dilatation in circumferential direction results only from remodeling (i.e., the dynamics of F i r governed by (17)) and can thus easily be overruled by sufficient growth (i.e., wall thickening) driven by F i g if the gain parameters k co σ and k sm σ are chosen high enough. Therefore, anisotropic growth in thickness direction can easily restabilize a blood vessel after a small perturbation. The situation is completely different in case of isotropic growth. There, by definition, the growth deformation gradient F i g entails an equal dilatation in wall thickness direction and in circumferential direction. While the wall thickening decreases the average wall stress, the dilatation of the circumference increases wall stress by the same factor, recalling Laplace' formula. Noting that the dilatation in circumferential direction is, however, additionally promoted by the remodeling deformation gradient F i r , it can easily overrule the effect of wall thickening, entailing a continued (and potentially unstable) dilatation of the vessel. It is worth noting that also in case of isotropic growth one may eventually observe a restabilization because constituents not subject to growth and remodeling (like elastin) will take over a larger and larger fraction of the wall stress as the vessel keeps dilating, which may enable collagen and smooth muscle to reattain at some point a homeostatic stress state. However, our computational studies suggest that for isotropic growth such a restabilization may occur only for large gain parameters (that is a very high collagen production) and after a considerable dilatation of the blood vessel. That is, isotropic growth appears to make blood vessels at least much more susceptible to unstable dilatation than anisotropic growth in thickness direction. If one assumes that growth and remodeling should ensure stability and functional adaptation of blood vessels as efficiently as possible, this raises the question whether or why blood vessels should rely on isotropic instead of anisotropic growth. In fact, Matsumoto et al. report experimental observations (in particular in Table 2 of Matsumoto and Hayashi (1996)) which suggest that growth and remodeling of blood vessels in response to hypertension are not isotropic. Rather the ratio of wall thickness and vessel radius increases significantly, which supports the idea that anisotropic growth in wall thickness direction rather than isotropic growth may be an appropriate model for vascular growth and remodeling. Further studies are, of course, required to definitely confirm or reject this hypothesis.

Conclusions
Mathematical modeling of soft tissue growth and remodeling has been based so far mainly on two major approaches. The kinematic growth theory of Rodriguez et al. (1994) relies on a simple multiplicative split of the deformation gradient into an inelastic growth part accounting for local changes of mass and volume and an elastic part ensuring mechanical equilibrium and geometric compatibility. While conceptually simple and computationally efficient, this approach suffers, from an inherent inability to account for the separate growth and remodeling of several distinct structurally significant constituents. To overcome this limitation, Humphrey and Rajagopal introduced the so-called constrained mixture models (Humphrey and Rajagopal 2002) where growth and remodeling are conceptualized as a continuous degradation and deposition of differential mass increments. These constrained mixture models realistically mimic the situation in vivo where cells such as fibroblasts and smooth muscle cells are continuously secreting new extracellular matrix and degrading extant matrix at the same time, for example, by matrix metalloproteinases (MMPs). However, the implementation of classical constrained mixture models based on multi-network theory is involved and their computational cost is significantly higher than the one of the simple kinematic growth models. To combine the advantages of simple kinematic growth models and constrained mixture models, recently a new class of models was introduced, the socalled homogenized constrained mixture models (Cyron et al. 2016b). They share the realistic micromechanical foundation of the classical constrained mixture models and produce results that are in many cases nearly identical or at least very similar to the ones of classical constrained mixture models. Yet, their implementation is nearly as simple as the one of the kinematic growth models and their computational cost is similarly low. By Cyron et al. (2016b) introduced homogenized constrained mixture models with a main focus on theoretical aspects.
In this paper, we pointed out in detail how homogenized constrained mixture models can be used for anisotropic volumetric growth and remodeling in three dimensions. Using the example of a circular-cylindrical model aorta, we demonstrated that homogenized constrained mixture models can reproduce realistically both pathological growth (as observed in aneurysms) and adaptive growth in healthy vessels (e.g., in response to hypertension). A major advantage of the homogenized constrained mixture models developed in this paper is the natural incorporation of anisotropic volumetric growth. Constrained mixture models for volumetric soft tissue growth and remodeling published previously were limited to isotropic volumetric growth. However, experimental observations, for example, by Matsumoto and Hayashi (1996), rather suggest that arterial growth is anisotropic and is likely to occur (at least mainly) in thickness direction of the vessel wall. As discussed in Sect. 4.3, the reason may be that anisotropic growth can ensure the stability of blood vessels under perturbations more efficiently than isotropic growth. Therefore, the ability of the homogenized constrained mixture model developed herein to incorporate straightforwardly anisotropic growth may make it an ideal tool to study growth and remodeling in the vasculature. In this paper, we used only a particularly simple constitutive model for the arterial wall with the intention to keep constitutive modeling as similar as possible to previously published twodimensional studies and enable thereby instructive comparisons. Combining the framework for growth and remodeling developed herein with more sophisticated and realistic threedimensional constitutive models such as the ones used by Eriksson (2014 will be a natural next step.

Appendix 1
In mechanobiological equilibrium, mechanical equilibrium is satisfied and at the same time no growth and remodeling occur because the stress of each constituent subject to growth and remodeling equals the homeostatic value . In two-dimensional membrane models of blood vessels, it is easy to define for a given vascular geometry an initial state in mechanobiological equilibrium. The reason is that in the balance of momentum only the membrane stress rather than the Cauchy stress appears (cf. Equation (2.2) in ). So, one can first solve this equation to ensure mechanical equilibrium. Subsequently, one can separately vary the remaining parameters (in particular wall thickness and mass fractions) to ensure that the Cauchy stress of each constituent equals the homeostatic value. By contrast, in three dimensions the definition of an initial state in mechanobiological equilibrium is not straightforward. In Sect. 4, we applied, inspired by Gee et al. (2010), the following procedure to initialize the simulations in mechanobiological equilibrium: 1. First, the homeostatic stress (or equivalently stretch) of collagen and smooth muscle in fiber direction has to be defined. Herein, we adopt the values used already by Wilson et al. (2013) and presented in Table 1. The elastic stretches λ co e (t = 0) and λ sm e (t = 0) of collagen and smooth muscle in the initial configuration are then computed so that they correspond to these homeostatic stresses (or are equal to the chosen homeostatic stretches). 2. Second, the elastic in vivo stretches of elastin in axial and circumferential direction have to be defined. Herein, we adopt the values used already by Wilson et al. (2013) and presented in Table 1. These stretches are equal to the initial stretch of elastin in axial direction λ el(90 • ) e (t = 0) and circumferential direction λ el(0 • ) e (t = 0). 3. Third, the elastic stretch λ el⊥ e (t = 0, r ) of elastin in wall thickness direction in the initial configuration at time t = 0 is computed. The stretch depends on the distance r from the cylinder axis and is computed such that the Cauchy stress σ ⊥ (r ) of the whole constrained mixture in wall thickness direction linearly increases from the value σ ⊥ (r = R) = −p at the inner radius of the cylinder to the σ ⊥ (r = R + H ) = 0 at the outer radius of the cylinder. That is, we use the condition With given strain energies from Sect. 3.3, given initial mass fractions (cf . Table 1), given initial (homeostatic) stretches of collagen and smooth muscle (cf. 1)) as well as given initial axial and circumferential stretch of elastin (cf. 2)), σ ⊥ in (40) becomes a function of only the parameter λ el⊥ e (t = 0, r ), which can thus be computed from (40). Note that the only other so far unknown parameter μ el 2D does not appear in (40) because it affects only the axial and circumferential wall stress but not the radial one due to the two-dimensional elasticity of (32). Equation (40) is solved for λ el⊥ e (t = 0, r ) by applying a Newton-Raphson method. 4. Fourth, the material parameter μ el 2D (r ) is defined such that it ensures a constant circumferential Cauchy stress of elastin despite the variation of λ el⊥ e (t = 0, r ) in radial direction. This ensures mechanical equilibrium for the λ el⊥ e (t = 0, r ) computed according to 3) and requires with σ i(0 • ) (t = 0) the Cauchy stress of the ith constituent in circumferential direction in the initial configuration at time t = 0. Equation (41) has to hold for any r ∈ [R; R + H ], and it is solved analytically for the material parameter μ el 2D (r ). Note that all parameters that are required to compute the Cauchy stresses σ i(0 • ) from the strain energies defined in Sect. 3.3 are known at this point except for μ el 2D . Also note that, as the elastic radial prestretch of elastin λ el⊥ e (t = 0, r ) varies in wall thickness direction, also the material parameter μ el 2D (r ) has to.
An initially stress-free thick-walled homogeneous cylinder on which an internal pressure is imposed exhibits in general a parabolic stress profile both of the circumferential stress σ (0 • ) and radial stress σ ⊥ as illustrated in Fig. 7a. By contrast, the prestressing procedure described in this appendix leads to a configuration with constant circumferential stress σ (0 • ) and linearly increasing radial stress σ ⊥ (cf. Fig. 7b) in case of uniform mass fractions and fiber orientations throughout the wall. The real stress profile over the wall thickness in arteries is currently not exactly known so that the prestressing procedure used herein should be understood as a simple way to initialize our simulations in a state of mechanobiological equilibrium rather than as a physiologically accurate procedure.

Appendix 2
In this appendix, we provide the first and second derivatives of the strain energy functions W co , W sm , and W el in (27), , and (36) with respect to the total Green-Lagrange strain E. For each constituent, we assume a multiplicative split of the deformation gradient into an elastic part F e and an inelas- Circumferential (σ (0 • ) ) and radial (σ ⊥ ) stress in a thick-walled cylinder from the inner radius R to the outer radius R + H in an initially stress-free cylinder under internal pressure p (left) compared to the prestressed configuration produced by the procedure described in this appendix (right) tic part F gr (cf. (13)) and the partial derivative with respect to the total strain is computed assuming that F gr is constant so that a variation of E translates into a variation of F e only. Note that in this appendix we omit superscripts in F gr and F e , given that in each equation it is evident to which constituent in the constrained mixture they refer. For collagen in (27), we have ∂ W co ∂E = 2 k co 1 (I a − 1) e k co 2 (I a −1) 2 F gr a 0 2 a 0 ⊗ a 0 where a 0 is the unit direction vector in reference configuration of the collagen fiber family. The second derivative of the strain energy function is ∂ 2 W co ∂E∂E = 4k co 1 1 + 2k co 2 (I a − 1) 2 e k co 2 (I a −1) 2 F gr a 0 4 a 0 ⊗ a 0 ⊗ a 0 ⊗ a 0 .
The stress and elasticity of smooth muscle fiber families aligned with the unit vector a 0 in reference configuration are governed by ∂ W sm pas ∂E = 2 k sm 1 (I a − 1) e k sm 2 (I a −1) 2 F gr a 0 2 a 0 ⊗ a 0 , ∂ 2 W sm pas ∂E∂E = 4k sm 1 1 + 2k sm 2 (I a − 1) 2 e k sm 2 (I a −1) 2 F gr a 0 4 a 0 ⊗ a 0 ⊗ a 0 ⊗ a 0 , For elastin we have, according to (36) where C e is the elastic Cauchy-Green deformation tensor of elastin, F gr its inelastic deformation gradient, A 0 = F −1 A gr F −T and the special tensor product is defined such that in index notation (A B) i jkl = A ik B jl + A jk B il /2. The three-dimensional neo-Hookean contribution of elastin leads to and with C gr = F T gr F gr and the standard Cauchy-Green deformation tensor C = F T F. The volumetric penalty term ensuring an isochoric elastic deformation of elastin leads to ∂ W el vol ∂E = 2ε |C e | 1/2 |C e | 1/2 − 1 C −1 and ∂ 2 W el vol ∂E∂E = 4ε |C e | 1/2 |C e | 1/2 − 1 2 C −1 ⊗ C −1 − |C e | 1/2 − 1 C −1 C −1 .
The above equations (42), (44), (46), (48), (50), and (52) can directly be used in (9) to compute for given reference mass densities of the different constituents the second Piola-Kirchhoff stress at each point. The elasticity tensor of the constrained mixture that is required to compute the tangent stiffness matrix is, as usual, given by where the second derivatives of the strain energies W i can be computed using (43), (45), (47), (49), (51), and (53).

Appendix 3
Time integration of growth and remodeling in a finite element scheme requires updating in each time step both the inelastic deformation gradients F i gr at each Gauss point and the total displacements at each node. For our homogenized constrained mixture model, we tested two time integration schemes. In both we discretized time by a series of n discrete points in time t = t 1 , . . . , t n with t 1 = 0. The distance between two subsequent points in time was t = t k+1 − t k . The total deformation gradient and inelastic deformation gradients at time t k are F t k and F i gr t k , respectively. In the first, explicit time integration scheme we compute at time point t k from a given total deformation gradient F t k and inelastic deformation gradient F i gr t k via the evolution equations (16)- (18) and (25) the rateḞ i gr and approximate F i gr t k+1 = F i gr t k +Ḟ i gr F i gr t k , F t k t. Then we compute from this F i gr t k+1 , using the balance of linear momentum (4) and the equations provided in Appendix 2, iteratively the mechanical equilibrium configuration at t k+1 , that is, F t k+1 .
In the second, implicit time integration scheme we first compute an estimate of F i gr t k+1 by solving (iteratively) the implicit equation F i gr t k+1 = F i gr t k + F i gr F i gr t k+1 , F t k t, whereḞ i gr is again computed using the evolution equations (16)-(18) and (25). Subsequently, we compute an estimate of the mechanical equilibrium configuration at time t k+1 , that is, F t k+1 , by solving (iteratively) the balance of linear momentum (4), using the equations provided in Appendix 2. Then we com-pute an updated estimate of F i gr t k+1 based on our estimate of F t k+1 instead of F t k . This updated estimate of F i gr t k+1 is used to update also the estimate of F t k+1 , and these updates of F i gr t k+1 and F t k+1 are iterated until F has converged.
Both in the explicit and implicit time integration scheme, the balance of linear momentum (4) has to be solved iteratively to compute F t k+1 . To this end, Newton-Raphson iterations can be used. Note that the tangent stiffness matrix used in these iterations can be directly based on (54) in Appendix 2 in an explicit time integration scheme. In an implicit time integration scheme, however, to obtain quadratic convergence of the Newton-Raphson iterations, one has to add to the second derivatives of the strain energies in (54) an additional fourth-order tensor accounting for expected changes of the inelastic deformation gradient through the iterations of the displacement field. This additional fourth-order tensor can be computed for the ith constituent as The computational results shown in Sect. 4 were all computed using the above-described implicit time integration scheme. Note that explicit and implicit refer in the above discussion to the update of the inelastic deformation gradient only. The balance of linear momentum is always obtained via the solution of an implicit system of equations.