Volumetric growth of soft tissues evaluated in the current configuration

The growth and remodelling of soft tissues plays a significant role in many physiological applications, particularly in understanding and managing many diseases. A commonly used approach for soft tissue growth and remodelling is volumetric growth theory, introduced in the framework of finite elasticity. In such an approach, the total deformation gradient tensor is decomposed so that the elastic and growth tensors can be studied separately. A critical element in this approach is to determine the growth tensor and its evolution with time. Most existing volumetric growth theories define the growth tensor in the reference (natural) configuration, which does not reflect the continuous adaptation processes of soft tissues under the current configuration. In a few studies where growth from a loaded configuration was considered, simplifying assumptions, such as compatible deformation or geometric symmetries, were introduced. In this work, we propose a new volumetric growth law that depends on fields evaluated in the current configuration, which is residually stressed and loaded, without any geometrical restrictions. We illustrate our idea using a simplified left ventricle model, which admits inhomogeneous growth in the current configuration. We compare the residual stress distribution of our approach with the traditional volumetric growth theory, that assumes growth occurring from the natural reference configuration. We show that the proposed framework leads to qualitative agreements with experimental measurements. Furthermore, using a cylindrical model, we find an incompatibility index that explains the differences between the two approaches in more depth. We also demonstrate that results from both approaches reach the same steady solution published previously at the limit of a saturated growth. Although we used a left ventricle model as an example, our theory is applicable in modelling the volumetric growth of general soft tissues.


Introduction
The interactions between living organs and the bio-environment play essential roles in regulating pathological or physiological growth. It has been experimentally demonstrated that environmental factors, such as the chemical, mechanical or genetic stimulus, could induce growth and remodelling (G&R) processes in living organs. Living organs can re-shape themselves, reset their constituents' growth (or turnover) rates, and develop volumetric and mass changes to adapt to pathological or physiological changes in the bioenvironment. G&R has been observed in different forms. For instance, nutrition concentration will be increased around the tumours, which helps to accelerate cell division and lead to the local growth of tumours at the tissue level (chemical factors). Embryonic or young organs can continuously restructure themselves to develop required functions. Mature organs are expected to stay in a relatively stable living state and serve as fully functional; however, pathologically, the dynamic impact of bio-environmental changes will induce a quick remodelling process to renovate the functional tissue. For instance, an embryo heart can continuously and instinctively develop its heart structure due to genetic factors, and tumours may grow independently of mechanical factors (Volokh 2006). However, a mature heart seems to grow along the direction of principal tensile stress due to mechanical factors (Taber and Eggers 1996). There are two typical examples of the maladaptive G&R in the heart. One is known as eccentric hypertrophy in response to chronic * X Y Luo xiaoyu.luo@glasgow.ac.uk 1 School of Mathematics and Statistics, University of Glasgow, Glasgow, UK volume overload, in which increased diastole wall stress leads to the addition of sarcomeres in series, associated with ventricular dilation. The other is eccentric hypertrophy due to pressure overload, whereby an increased systolic wall stress leads to the addition of sarcomeres in parallel, resulting in wall-thickening. Despite the importance of G&R in disease progression, the principles governing those mechanisms are not fully understood. Improving this field's knowledge will have considerable potential for optimizing clinical treatments to save more lives efficiently. However, although growth has been of research interests for a few decades, and many different G&R theories exist but these mostly fall into two categories. The first one is based on variable material points, and the second one is based on fixed material points. The most well-known approach in the first category is known as the "constrained mixture" theory, initially developed by Humphrey and Rajagopal (2002), in which living tissues are considered to be composed of different constituents which continuously turnover with mass changes during the life cycle. They postulate that each constituent has a "preferred state" and that the reference (natural) configurations for each constituent exist and will be whatever allows this preferred state to be achieved. This framework considers the reproduction and removal of each constituent, so the reference configurations can be different even for the same constituents. In essence, one can add or remove material points to change the reference configurations. Many researchers, e.g. Baek et al. (2006); Alford et al. (2008); Valentin and Humphrey (2009), and Watton et al. (2004); Watton and Hill (2009) developed their models in this category to study the G&R of arteries and aneurysms. The difficulty of using this theory is the lack of experimental data for historical changes of all the constituents during the tissue growth and it is almost impossible to track different natural configurations from experiments. Stimuli from combined effects of various sources, including mechanical, thermal, electrical, and genetic information, may simultaneously contribute to the successive growth processes in living organs over different time scales, raising the complexity of identifying the contributions of individual environmental factors.
The volumetric growth theory is a typical example of the second category, following the idea of the kinematics of multiplicative plasticity (Bilby et al. 1955;Kerckhoffs et al. 2012). In this theory, it is assumed that there exists one reference (natural) configuration for all the constituents, and any incompatible growth due to growth of different constituents can be absorbed by the residual stresses. Taber (1998) and Rachev et al. (1998),among others, have employed this concept to model arteries and their response to hypertension in maturity. Following a growth law introduced by Taber (1998); Kerckhoffs (2012) used a strain-driven law to model post-natal cardiac growth, in which a cumulative growth is considered using multiplicative decomposition of a consecutive sequence of growth and elastic deformation tensors. Moreover, Kerckhoffs et al. (2012) modified the growth law to exclude the "unbounded" growth with a growth multiplier limiting function. Göktepe et al. (2010a) also developed a framework for using the bounded stress-and strain-driven growth laws to study pathological and physiological growth processes in the living heart. One of the debating points in this approach is whether a growth field is compatible. In addition, since soft tissues are essentially fibre-reinforced materials, growth theories that describe changes in the fibre structure have also been put forward, either by directly tracking the micro-fibre structure remodelling (Driessen et al. 2003) or indirectly by including the density growth where the volume of the structure remains constant but its density varies (Eriksson et al. 2014). A few studies considered the evolution of the fibre distribution as well during G&R (Kroon et al. 2007;Rouillard and Holmes 2012;Zhuan et al. 2019).
Of the two categories, the volumetric growth theory has been widely used due to its simplicity. However, a longstanding issue of the volumetric growth theory is that most existing volumetric growth studies assumed that growth occurs in the natural (reference) configuration (Budday et al. 2014;Göktepe et al. 2010a). Under this assumption, the volumetric growth is modelled by introducing a socalled growth tensor defined from the reference configuration (Rodriguez et al. 1994;Hsu 1968;Kerckhoffs et al. 2012;Kerckhoffs 2012;Klepach et al. 2012), so that the overall deformation tensor can be decomposed as where g is the pure growth tensor, and e is the pure elastic deformation tensor. For fibre-reinforced soft tissue modelling, if the fibre structure in the reference configuration is represented by the fibre, sheet, and normal directions 0 , 0 , 0 , then a commonly used growth tensor can be expressed as  where j ( , t) , (j = f , s, n) , are called the growth multipliers. Note that growth is assumed to be along the principal cylindrical directions, and are functions of the position vector of the reference configuration and time t. These are smaller or greater than one when the elastic body shrinks or grows with respect to the reference configuration.
If growth is compatible, then no residual stress is generated. This process is discussed in sect. 2.1.1. The evolution of the growth (multipliers) can be either considered as stress-driven (Taber 1998;Lubarda and Hoger 2002;Taber and Chabert 2002;Hariton et al. 2007), strain-driven Driessen et al. 2004), or combined stress-and strain-driven (Taber and Eggers 1996). One important aspect of G&R theories is their thermodynamics consistence, which has been discussed in details by Epstein and Maugin (2000). For growth from the reference configuration, it appears that the Mandel stress plays an important role in thermodynamically consistent formulations of inelasticity. It is for this reason that Mandel stress is often used in stress-driven growth laws. Amar and Goriely (2005) and Goriely and Amar (2007) moved a step further and studied growth from the current configuration without releasing the residual stress, albeit for simplified cases, such as symmetrical spheres (Goriely and Amar 2007). However, living organs cannot always be so simplified. More recently, Genet et al. (2015) introduced a growth involution law in the updated stress-free configuration. However, they still used the growth tensor defined in the fixed reference configuration, and the strain energy function remained the same during growth. Therefore, it is unclear if their approach is thermodynamically consistent.
In reality, soft tissue G&R is a continuous process, so growth must occur in the current (loaded and residually stressed) configuration. In other words, given the geometrical constraints, pure growth generally induces material incompatibility (Skalak et al. 1996). However, a theoretical framework for an incompatible growth in the current configuration is missing, although some efforts have been made in this direction by studying growth evolution from an updated reference configuration after each incremental growth. For example, an inhomogeneous volumetric growth was studied using a three-dimensional simulation of the heart (Kroon et al. 2009). However, in these studies, it was assumed that residual stress is released after each incremental growth, and that further growth always starts from the updated, yet stress-free configurations. The reason for people to adopt this assumption is that the updated natural configurations can be determined directly from the cumulative incremental growth (or growth history) which leads to a similar computational method to the growth models from the reference configuration (Budday et al. 2014;Göktepe et al. 2010a). However, the cost of this convenience is to force growth to be compatible at all times, which is a major model limitation. Amar and Goriely (2005) and Goriely and Amar (2007) moved a step further and studied growth from the current configuration, albeit for simplified cases, such as symmetrical spheres.
In particular, Goriely and Amar (2007) employed an incremental theory to propose a continuous description of the sequential growth from the current configuration. Although they started with general problems, progress was only made when the growth tensor and the deformation gradient tensor are both assumed diagonal.
More recently, Genet et al. (2015) introduced a growth evolution law in an updated stress-free configuration. However, they still used the growth tensor defined in the fixed reference configuration, and the strain energy function remained the same during the growth. Therefore, it is unclear if their approach is thermodynamically consistent.
In this paper, we propose a new framework of volumetric growth evaluated in the current configuration, which is thermodynamically consistent, and without releasing the residual stress or imposing geometrical and deformational restrictions. In addition, we shall show that our theory can recover the results of Goriely and Amar (2007) under the same simplifications they used. We also establish the relationships between the strain energy functions defined in different configurations, i.e. in the natural, current (stressed), and the evolving, grown, but stress-free configurations. This is different to previous work done by groups, e.g. (Ogden and Saccomandi 2007;Holzapfel and Ogden 2010;Gower et al. 2017;Agosti et al. 2018), who derived the constitutive laws from current stressed configurations for soft tissue mechanics, but not applied to G&R processes.
We illustrate our idea using both a left ventricle (LV) and a cylinder model, which admit an inhomogeneous growth in a residually stressed current configuration. We also show that the total cumulative growth from the reference configuration, the expression of which is prescribed in previous volumetric growth theories, can be derived based on this new framework. This cumulative growth is not affected by elastic stretches but is dependent on elastic rotations. In other words, the cumulative growth tensor is a function of loading history, and hence, cannot be postulated a priori. The strain energy function with respect to the updated reference configurations also changes with growth. In this new approach, the stress-driven growth law is taken to be a function of the Cauchy stress tensor (also known as the true stress). Finally, we compare the residual stress of the LV model using growth laws defined in the reference and current configurations, respectively, with published experimental measurements (Costa et al. 1997;Omens et al. 1993), and show that our approach leads to qualitative agreements with the measurements.

Volumetric growth theories
We take the accepted hypothesis that pure growth does not induce any elastic deformation [H1]. In other words, pure growth and elastic stretching of a body are two independent events. We also assume that the strain energy of a body at any time is a function of the total deformation gradient, which can be a combination of elastic and residual-stress induced deformations [H2] (Ogden 1997).

Growth from the reference (natural) configuration
We first introduce the concept and notations of existing volumetric growth theories based on the reference (natural) configuration.

Compatible growth
Let and be the position vectors of a material point in the reference and current configurations, B 0 and B t , respectively. The idea of growth from the reference configuration is illustrated in Fig. 1, where a pure (compatible) growth g takes an elastic body from B 0 to an (imaginary) grown and stress-free configuration B g as g = 0 0 , assuming that the map 0 ( 0 ) exists and is differentiable. Then, after loading the body deforms into the current configuration B t with elastic deformation gradient e = 0 .
The overall deformation gradient is Following (Ogden 2003), for incompressible materials that undergo volumetric growth, we must have the connection where W 0 is the free energy per unit volume in B 0 , W g is the elastic strain energy function per unit volume in B g , J A = det , and J F e = det e . Since W 0 ( g ) ≡ 0 due to the hypothesis [H1], we simply have where we have used J A = J g = det g , and J F e = 1.
(3) = e g . (4) The Cauchy stress in B t can be determined either through the total deformation from B 0 , or through the elastic deformation e from B g , Equation (6) can be proved by making use of (5) and (3), i.e.

Incompatible growth
We use the pure growth to describe the process so that if the growth tensor used to designate the growth which would be locally observed at a point of the tissue with the small region around it could be grown in isolation, then this leads to an incompatible strain field, and residual stress is required to keep the tissue intact (Skalak et al. 1996). This process is shown in Fig. 2, where the grown and stress-free configuration B g is made compatible via an elastic deformation into B with residually stress, before the loading-induced e maps B into B t . As shown in Fig. 2 and [H2], this process can also be viewed as a growth-induced residual deformation from B g to B followed by e .
We now focus on incompressible materials and choose B g as the new (evolving) reference configuration, with a strain energy function W g . The total elastic deformation gradient from B g to B t is (7)

Fig. 1
Compatible growth from the reference configuration, B 0 . A pure and compatible growth g takes the elastic body from B 0 into B g , and is then deformed via the loading-induced elastic deformation e into the current configuration B t . The total deformation gradient from B 0 to B t is denoted as Fig. 2 A pure growth g from B 0 takes the body to an incompatible and stress-free state in configuration B g . The deformation assembles the body into a compatible, but residually stressed stated in configuration B . Loading induced deformation e then maps it onto the current configuration B t . The process can be equally viewed as −1 g from B g to B 0 followed by from B 0 to B t , or from B g to B t directly via the elastic deformation E As before, the Cauchy stress in B t can be determined either from the reference configuration B g , where W g is the elastic strain energy function related to the configuration B g .
or from B , where the strain energy function W from B is a function of both e and residual stress , and p, p are the Lagrange multipliers to ensure incompressibility of the material. The identity ( p = p ) has been proved by (Ogden 2003). Henceforth, we only use p.
The equivalence of (9) and (10) can be similarly proved using the connection, Substituting (5) and (11) into (10), and making use of [H1], so that ∕ e = , we obtain (9). When evaluated in B with e = and E = , the Cauchy stress simply becomes the residual stress , i.e.
The residual stress and the Cauchy stress must obey their corresponding equilibrium equations, The boundary conditions are zero-loading for and loading for , respectively. This enables us to solve for and e , and hence, and .

Growth evaluated in the current (stressed) configuration
We now propose a new theory that enables growth law to be evaluated in the current (stressed) configuration. We denote such a configuration as B 1 . As shown in Fig. 3, B 1 is reached after the body is deformed from the reference configuration B 0 .
Let 1 =̄ e 1 indicate this elastic deformation. Since the line elements of material transform according to d = 1 d , we have the incremental connection, (13) div = or div = . (14) If the stressed configuration B 1 is now chosen as the reference configuration then the right-hand side of (14) becomes (d ) , where is the value of 1 in this configuration. Here, 0 is the unit tensor before the initial configuration B 0 : 0 = I . Since d and hence (d ) is independent of the reference configuration, we obtain the connection (Ogden 1997) Now consider an incremental deformation caused by the pure growth ̄ g 1 and ̄ 1 , which take the body into the incompatible configuration B g 1 first and then assemble it into a compatible configuration B 1 . If the growth displacement is "small" for each in B 0 so that terms of order ≥ | | 2 are negligible in comparison with those of order | | then we refer to 1 =̄ 1̄ g 1 as an incremental (growth) deformation from B 1 to B t . This approach is the classic growth incremental procedure (i.e. in Amar and Goriely (2005); Göktepe et al. (2010a)). Goriely and Amar (2007) also employed the similar approach but included the second-order deformation increments. Since can also be written as As in Sect. 2.1.2, the total elastic deformation gradient E 1 can be equally achieved either along the top path , which can be derived from ̄ g 1 , maps the reference configuration B 0 into the grown and incompatible (stress free) configuration B g 1 . B g 1 is made compatible by 1 and becomes B 1 with the residual stress 1 . From B 1 the body is deformed to B t by the elastic deformation e 1 . The process . Then, the body undergoes a small incremental The process can be equally viewed as −1 We remark that for growth from the current configuration, it is not easy to follow the top path in Fig. 3, since we do not introduce the strain energy function with respect to B g 1 or remove the external loading. On the other hand, if we follow the lower two paths, we have the situation already discussed in Sect. 2.1.2. In other words, we always compute from (18), and not from (16). However, to use (18), we must first establish the connection between ̄ g 1 and g 1 . Since pure growth is independent of the elastic stretch (H1), the only difference between ̄ g 1 and g 1 is due to the rotation of ̄ Since ̄ e 1 can be easily computed, given ̄ g 1 , g 1 is known. The rest of the analysis simply follows the approach of Sect. 2.1.2. The Cauchy stress in B t is The strain energy function W 1 with respect to B 1 satisfies The residual stress 1 is simply the Cauchy stress evaluated in B 1 when the external loading is removed, i.e. e 1 = , Again, 1 and 1 both obey the corresponding equilibrium equations, and the appropriate boundary conditions. From these, we solve for 1 and e 1 , and 1 and 1 .

Subsequent growth from stressed and grown configuration
Subsequent growth can occur for k steps of growth and deformation in B k , as shown in Fig. 4. In general, external loading may also change after every step. Now with an new incremental growth k from B k , the total deformation gradient from B 0 to B t becomes Here, G k is the total cumulative pure growth from B 0 to B g k , computed from the previous cumulative pure growth Again, to follow the path B g k -B k -B t , we must first find the connection between ̄ g k and g k , so that (26) can be evaluated. The deformation gradient between B k and B gk is , then we have the connection Fig. 4 General roadmap after k steps of continuous growth and deformation. The deformation gradient from B 0 to B t is , the incremental growth from the loaded configuration B k is k =̄ k̄ g k . The total cumulative pure growth tensor B 0 to B g k is G k . B g k is the grown, stress free but incompatible configuration, and B k is the compatible and residually stressed configuration where r k is the rotation of r k .
The procedure is now the same as before. So we give the general expressions. The Cauchy stress in B t is where Similar to (22), the residual stress k is Again, k and k both obey the corresponding equilibrium equations and the appropriate boundary conditions. From these, we can solve for k and e k . At this point, all we need is an incremental growth law for ̄ g k .

Incremental growth law
If we consider continuous growth and deformation, then we can write g k = g t , similarly, G k = G t . For simplicity, we assume that the growth tensor is diagonal for a fibrereinforced material. This assumption is widely accepted by the community and seems to hold for the left ventricle (Li et al. 2021). Consequently, all the incremental growth tensors are also diagonal. Note it is not essential to assume a diagonal incremental growth tensor in this framework. A more general growth law including off-diagonal elements can equally work if one has experimental data to support the expression. However, the assumption is used to illustrate the framework more clearly, and to compare our results with those from previous work. Let ̄ g t = d be the incremental growth tensor, we write where ( t , t , t ) is the fibre structure of the material at time t, which is related to the fibre structure at t = 0 via and ̇f ,̇s,̇n in (31) are the rates of change of the growth multipliers in the t , t , t directions, respectively, which can be written as ̇ ( , t) = ̇f ,̇s,̇n . Note = diag( G t ) , so represents the total cumulative growth vector of the elastic body with respect to the reference configuration B 0 .
We now assume that growth can be either stress or strain driven. Using a similar evolution law as in (Lubarda and Hoger 2002;Göktepe et al. 2010a), we write where t and t = 1 2 ( − −T −1 ) are the Cauchy stress and Eulerian strain tensors in B t , respectively.
The growth threshold is a scalar function of either the Cauchy stress tensor or the Eulerian strain tensor, depending if growth is stress or strain driven. In other words, growth is activated once the stress or strain during physical activity exceeds the physiological threshold level. If is negative, then there is no growth. l= l f , l s , l n is a limiting or scaling vector function since growth cannot continue infinitely. In this paper, we follow Lubarda and Hoger (2002); Göktepe et al. (2010a) to choose where jmax are the maximum values of the growth multipliers j , and j , j are growth parameters. The choice of and l j ( ) in (33) will ensure the growth rate changes smoothly until the growth multiplier j has reached its maximum value jmax . The nonlinearity and the speed of the growth are governed by j and j , respectively.

The strain energy function
We choose the grown and stress-free configuration B g t as the reference configuration which evolves with time. We also adopt a modified Holzapfel-Ogden (HO) model for the fibre-reinforced material (Holzapfel and Ogden 2009) where I 1 , I 4 are invariants of E t = T E t E t , the right Cauchy-Green tensor with respect to B g t , E t is the total elastic deformation from B g t to B t , g t and g t are angles between g t and 0 , and g t and 0 in the 0 − 0 plane, respectively, and the fibre direction in B g t changes according to And a, b a 1 , b 1 are material parameters.

Finite element model and numerical algorithm
The LV diastolic dynamics with growth from loaded configuration is simulated using the finite element package FEAP 8.3, with 2100 8-node hexahedral elements and 2375 nodes (Fig. 5). In addition to the equilibrium equations and boundary conditions, where we fix all the nodes on the base, we also need to solve the evolution equation (33). Since the continuous growth occurs at a time scale much slower than the time for the system to reach the equilibrium, the equations are quasi-time-dependent.
At the current time t, after k steps of growth. The growth evolution equation (33) predicts the incremental growth At every time step, we use a mixed variational approach in FEAP to solve the mechanical behaviour of the nearly incompressible soft tissue that is free from locking. For each element, FEAP uses the Hu-Washizu variational principle where ̂ , ̂ are the element ordered vectors of stress and strain tensors, u is the displacement vector, ℂ (2) is the ordered matrix of tangent moduli tensor (4) (see (Taylor 2014) for details of the definition of an ordered matrix), t is the traction, ∇ (s) is the symmetric part of ∇ , b is the traction on the force boundary , b is the displacement vector on the displacement boundary . The computational algorithm for the volumetric growth is given below. 5 a The LV geometry with 28mm long axis, an internal radius of 5mm and external radius of 10mm at the base, and a block cut from the LV wall. The ratio between the LV wall thickness and the internal radius is chosen to be one following Omens and Fung (1990). The basis vectors at the reference configuration are ( 0 , 0 , 0 ) for local coordinates, where 0 , 0 , and 0 are the local circumferential, longitudinal and transmural unit vectors. The basis vectors at the reference configuration are ( , , ) for global Cartesian coordinates, with origin O at the LV apex. b The fibre structure through the thickness of the LV wall. c Five longitudinal-circumferential sections through the wall thickness. Collagen fibres lie in the 0 -0 plane

Test examples
We now use the following examples to illustrate our approach. Henceforth, we use Growth-B 0 to refer to the volumetric growth from the reference configuration, and Growth-B t to indicate the growth from the current configuration.

Geometry of the LV model
Myocardium is considered to be a fibre-reinforced material composed of collagen fibres and myocytes. A local coordinate system with the circumferential, longitudinal, and transmural basis vectors ( 0 , 0 , 0 ), is introduced to describe the layered fibre structure within the ventricular wall, as shown in Fig. 5. An idealized half ellipsoid geometry is used for a rat LV. Global Cartesian coordinates (X, Y, Z) are used to describe material points in the undeformed reference configuration, with the corresponding basis vectors denoted { X , Y , Z } . Note that Let the myofibre architecture be described by a "fibre-sheet-normal" system ( 0 , 0 , 0 ) in the reference configuration, B 0 (Wang et al. 2013). Here, we assume that the fibre direction 0 always lies in the 0 -0 plane, the sheet direction is transmural, and the sheet-normal 0 = 0 × 0 .
We consider two cases of growth in the myocardium, one is the isotropic growth of athletic's heart, and the other is the growth following a pathological cardiac dilation. For motivation and physiology, please refer to Göktepe et al. (2010a), who first considered these cases. In both cases, to simulate the LV growth under loading, we use the strain energy function defined in (35), with parameters chosen to be a = 2.28 kPa, b = 1.8 . We also impose a constant internal pressure of 12 mmHg from t = 0 throughout the growth.
We note in this special case, the fibre structure ( t , t , t ) in the updated stress-free configuration B k is the same as the initial fibre structure ( 0 , 0 , 0 ).
For Growth-B 0 , we follow Göktepe et al. (2010a) by assuming the growth is stress driven, using the growth tensor defined in (2), together with the evolution law The corresponding experimental measurement of e ff from a canine heart mid-anterior wall by Costa et al. (1997) is also plotted in (a) for comparison where Here is the Mandel stress tensor, = , = T , is the second Piola-Kirchhoff stress tensor, and M crit is a critical stress value. The parameters are chosen as For Growth-B t , we have the growth tensor (40) and its evolution law defined in B t , and (33) 1 and (34), with and similarly choose with all other parameters being the same as in Growth-B 0 . Figure 6 shows that both models display ventricular dilation and wall thickening.
However, the growth patterns are clearly different. In Growth-B 0 , most of the growth is focused on the apex area, which is opposite to that of the Growth-B t model. Away from the apex area, the growth is also very different. To make more detailed comparisons, we define a transmural path 1 (see Fig. 6) at the same location for both models.
We then compare the Eulerian residual strain and stress components, e ff and ff , along the mean fibre direction, calculated from the residually stressed configuration B t as shown in Fig. 4. The transmural distributions of e ff and ff at weeks 0 to 5 are shown in Fig. 7 for both models. We can see that the trends of these distributions are opposite in the two models. For Growth-B 0 , e ff is positive from the endocardial surface, and decreases transmurally; but for Growth-B t , e ff is negative from the endocardial surface, and increases transmurally. The trend of Growth-B t is supported by the published experiments and previous studies, e.g. Costa et al. (1997); Omens et al. (1993); Wang et al. (2014). In particular, the experimental data by Costa et al. (1997) measured from the mid-anterior part of the canine LV are drawn in Fig. 7 to illustrate the qualitative dis/agreements with the two models.

LV case 2: Pathologic cardiac dilation
In the second example, we consider dilated cardiac growth, which is often a result of myocardial infarction (heart attack). In this case, the heart responds to a volume overload while attempts to maintain the cardiac output at a physiological level (Cheng et al. 2006). This cardiac dilation can be described as a strain-driven, transversely isotropic, and irreversible growth. Again we follow Göktepe et al. (2010a) and model the process assuming For Growth-B 0 , the growth tensor is defined in B 0 as The strain-driven involution law is Göktepe et al. (2010a), where l( ) is the same as in (41)  For Growth-B t , we have the incremental growth tensor defined in B t , We use the same strain-driven law as (43), except the function is function of where crit f is the criteria value of elastic stretch along t . In this case, the fibre structure in either B t or B g changes with time, so we need to compute the updated fibre structures using (32) and (36). Other than this, we keep the same parameters as in (45) with crit f replacing crit f 0 . The distribution of for both models is compared in Fig. 9. Again, similar patterns are seen in each model, with more growth in the apex area for Growth-B 0 , and more growth away from the apex area for Growth-B t . We also see that compared to the isotropic growth in the previous example, in both models, this strain-driven cardiac growth (eccentric dilation) shows a more significant increase in the cavity size, but with a less noticeable increase in the wall thickness.
The residual strain e ff and stress ff along the Path 1 of the LV are shown in Fig. 10 for both models. Again, the two models have opposite strain distributions. Note that we do not see wall thinning of the LV wall in either model, which is frequently observed in clinic cases, e.g. (Hankiewicz et al. 2008;Venco et al. 1987), or modelled in (Zhuan et al. 2019). However, here no local growth or infarction zone is considered, which is responsible for wall thinning. Similar simulation and s = g and s = n ≡ 1.
conclusion are seen in (Göktepe et al. 2010b;Lee et al. 2015;Klepach et al. 2012). Similar trends of the maximum residual loop stress in Fig. 8 are observed in this case.

Saturated growth in a multi-layer cylinder
To understand the differences between the two approaches in greater depth, here we consider a stress-driven growth in a cylindrical model subject to internal pressure, similar to that studied by (Goriely et al. 2008), except we use four layers here to enable the cross-layer growth. This model also allows us to compare to the special case considered by (Goriely et al. 2008), where a constant growth tensor (saturated homogenized growth) is used.
The model geometry in terms of the polar coordinates (R, , Z) is given by where L is the LV length, and R (i) j , R (o) j denote the inner and outer radii of each layer ( j = 1, 2, 3, 4 ), respectively.
In terms of cylindrical polar coordinates (r, , z), the geometry of the deformed configuration is given by where is the axial stretch. Similar to the incremental growth used by Goriely et al. (2008), we choose the incremental growth tensor in Growth-B t to be and the growth tensor for Growth-B 0 to be with a simplified stress-driven growth involution law as in , and (48) (1, , 1), where and max are chosen to be 1 and 1.75, which are arbitrarily set to reach the limit (saturated growth) sooner.
Notice for this simple model, the rotation tensor caused by elastic deformation is always identity; therefore, the cumulative growth tensor G t of Growth-B t is the same as the growth tensor g in Growth-B 0 . Note in the approach by Goriely et al. (2008), no stress-or strain-driven law is used; the cumulative growth tensor for the saturated growth is simply given as where ∞ is chosen to be max here.
Solving the momentum equilibrium equation with compatibility conditions as in (Zhuan and Luo 2020), we obtain solutions of the problem. Figure 11 shows the transmural distributions of the residual hoop strain e and stress for the 4-layer cylindrical model at different growth times, using Growth-B 0 and Growth-B t approaches. Figure 12 plots the cumulative growth factor and the time history of the maximum value of (the history of the residual loop strain has a similar trend, not shown).
Note that in both models, increases monotonically with time. However, the trend of stress/strain is different; it increases monotonically with time for Growth-B 0 , whereas for Growth-B t it reaches a local maximum first before decreasing. When growth is saturated, Growth-B 0 and Growth-B t both reach the exactly same solution as in (Goriely et al. 2008).
To be more general, we can assume the cylinder is made of many layers of equal thickness d in the reference configuration. We consider the cumulative diagonal growth tensors for all the layers After the growth, the jth layer is deformed with the middle perimeter 2 (R (i) j + d∕2) j , and the radius (R (i) j + d∕2) j , where R (i) j indicates the inner radius of the jth layer. Similarly, the middle radius of the next layer is ( indicates the inner radius of the (j + 1) th layer. Material continuity (compatibility) requires that (59) can be rearranged to read If we consider infinitely many layers, i.e. taking the limit of d → 0, then (60) is simply (1, j , 1), j = 1, 2, 3, ...(from inner to outer walls).
where ′ is the derivative of with respect to R. This is the material compatibility condition. If (61) holds, then the local volumetric growth is compatible and does not induce residual stress. If (61) is true everywhere, then the whole configuration is residually stress free. Hence, we can introduce the measure of incompatibility as If G > 0 , the growth will induce a smaller residual hoop stress in the jth layer compared with the (j + 1) th layer. If G > 0 across the whole wall, then the hoop stress distribution will be from negative (compression) to positive (stretching). Likewise, if G < 0 , the growth will induce a greater residual hoop stress in the jth layer compared with the (j + 1) th layer. If G < 0 across the whole wall, then the hoop stress distribution will be positive to negative. Since the cylindrical model is subject to internal pressure, the distribution of the Cauchy hoop stress is from positive to negative in the absence of growth. With the given growth (58), we found G > 0 across the whole wall initially in the Growth − B t model, leading to a positive maximum residual hoop stress towards the outer-most wall, as shown in Fig. 7b. This is opposite to the distribution of the Cauchy hoop stress, and tends to even out the overall stress distribution. In other words, the residual stress should have a beneficial effect under physiological loading. As this trend continues, growth, driven by the total stress, will have to slow down at some stage, and the growth mismatch, hence, the residual stress, will start to reduce. This occurs at t=1.5 in Fig. 12.
On the other hand, in the Growth − B 0 approach, we found G < 0 across the whole wall initially, leading to the negative maximum hoop stress towards the outer-most wall, as shown in Fig. 7a. This distribution of the residual loop stress holds the same trend with the loading induced Cauchy stress. Hence, the overall stress increases, which in turn enlarges the growth mismatch. The (negative) residual hoop stress becomes more negative with time.
Since growth is set to be limited by max , giving sufficient time, all the material points eventually reach the same saturated growth in both approaches, so we have which leads to G < 0 everywhere. In other words, for a saturated growth, we have the same negative maximum hoop stress using both approaches, as predicted by Goriely et al. (2008).

Discussion
In this paper, we have developed a new approach that allows the growth of soft tissues to occur in the current loaded and residually stressed configuration, which is different to the traditional approaches that need to define growth tensor from the reference (unloaded and zero-stressed) configuration. Since all functional activities and remodelling processes in living tissues continuously occur and respond to bi-environmental signals, it is natural to assume that growth is induced in the current configuration to adapt to the time-dependent changes. Although we also assumed that the growth tensor is diagonal in the examples, our theory also works for more general expressions of the growth tensor. The striking differences in the current approach, Growth-B t , and the previous approach, Growth-B 0 , are demonstrated in two test cases. Namely, isotropic growth and pathologic dilation of myocardium. The most significant finding is that the estimated residual strain distributions post G&R have opposite transmural distributions using the two approaches. And the result of the Growth-B t model is supported by experimental observations. The other differences we noticed are that in the Growth-B t model, the overall residual fibre stress or strain seems to reach a peak first (e.g. about week 3), and then decreases with time. This is different to the Growth-B 0 approach, where the residual stress or strain increase monotonically with time until the set growth limit is reached, as shown in Fig. 8.
These differences are explained in detail using a cylindrical model. Indeed, the cylindrical model has the same characteristics as the myocardium models, in terms of the behaviours of the two approaches. Further, with this model, we are able to introduce an incompatibility index, which can be used to judge the trend of the transmural distribution of residual stress. It is also easier to demonstrate that at sufficiently long time ( t > 5 in this case), both models converge to the same saturated growth limit (black curve), as the growth factor of all regions reaches the set limit of max , which is the same solution of the constant growth studied by Goriely et al. (2008) for > 1. Fig. 9 Contours of the growth factor in the Case 2 after 0, 1, 3, and 5 weeks growth. Above: Growth-B 0 . Below: Growth-B t . In both cases, the LV cavity size increases more obviously than the wall thickness (a) (b) Fig. 10 As in Fig. 7, but for the Example 2 For the saturated growth in the cylindrical model, we have the same negative maximum hoop stress using both approaches, as predicted by Goriely et al. (2008). In other words, in the saturated growth, the residual stress does not help the soft tissue cope with the external loading, which is against the experimental observations. However, we remark that this limiting steady state is a result of the prescribed growth law, which may never be achieved in physical situations within a realistic time scale. Not only because the residual hoop stress distribution of the limiting case is not seen in experiments, but geometrical and time constraints of a living tissue/organ may not allow a homogenized (saturated) growth to occur. Clearly, in the two myocardium cases, none has reached the saturated growth. We remark that the cumulative volumetric growth G t from B 0 essentially plays the same role as the conventionally used growth tensor 0 g in (2). In other words, once we know G t , all the classical theories of G&R follow. However, it is impossible to define an evolution law for G t prior to deformation, since elements of G t are continuously rotated due to elastic deformation, as shown in (26) and (27). On the other hand, one can easily prescribe an evolution law for the incremental growth tensor g t . This explains why the updated growth theory by Goriely and Amar (2007) works for symmetrical geometries with diagonal deformation gradient and growth tensors, since deformation induced rotations do not make any differences. In this special case, a growth law for G t may be defined prior to deformation.
To be more specific, we now show that under the assumption that both the deformation gradient and growth tensors are diagonal, the residual stress tensor arising from the incremental growth tensor ̄ gk used by Goriely and Amar (2007)(i.e. " 2 1 " in their Fig.2), is the same as in our approach. Note the rotation of ̄ gk is an identity tensor due to its diagonality (i.e. " 2 ′′ ≡ )) and " 1 " is essentially the mapping from B g k to B g k in Fig. 4, using our notation this is k In other words, for this special case, we obtain the same results as in Goriely and Amar (2007). However, when either the growth tensor or the elastic tensor is not diagonal, these two approaches are different. Hence, for general G&R processes, our residual stress tensor ̄ gk k̄ −1 gk cannot be simplified to k −1 G k−1 and is dependent on the past loading history.
G&R modelling is useful even if we do not know the history of the tissue growth. Often we encounter situations when all we know is that residual stress exists. A growth theory can be used to establishes the relationship between the residual strain/stress and the history of tissue growth. Hence, with a given growth tensor, we can estimate the residual strain (Zhuan and Luo 2020). On the other hand, by estimating the residual strain experimentally, e.g. via the opening angle methods, we can determine the cumulative growth tensor at the time when the residual stress is released. By considering growth in the current configuration, we also naturally included the effects of changed fibre orientations in the soft tissue. Indeed, the mean angle of the fibre distributions changes after continuous growth, which in turn updates the strain energy function.
Finally, we would like to remark that in this work a simplified HO model is used and the tuning factors in the test cases are not based on experiments. However, we believe these do not affect the qualitative behaviour of the outcome.
Although our focus is on the myocardium tissue, the developed G&R framework is applicable to all soft tissues, and could pave the way for more realistic growth modelling. In reality, the limit of tissue growth should be dictated by local maximum stress/strain or physiological homeostasis, given boundary conditions. It is seldom that a biological organ can reach a saturated global growth. We remark that finding a suitable growth limit based on experimental evidence is crucially important for practical applications. This research topic remains a challenge.

Conclusion
In this paper, we proposed a new theory of volumetric growth which depends on fields evaluated in the current configuration. In contrast to most previous models that have to define soft tissue growth in the natural configuration, we allow such a growth to be evaluated in the deformed and loaded configurations, without imposing any geometrical restrictions. We illustrated our idea using a simplified left ventricle model, that admits inhomogeneous growth in residually stressed and loaded configurations. We then compared our residual stress distribution with a typical previous volumetric growth model in which the growth tensor is defined in the natural configuration. Our results show that the new framework leads to a qualitative agreement with experiments, in contrast to the results from the previous model. In addition, using a multi-layer cylindrical model, we explained the nature of the differences in these two approaches using an incompatibility index, and demonstrated why, after an infinitely long time, both approaches attain the same homogenized/saturated growth state as identified by Goriely et al. (2008). We further observed that this limit may never be reached in soft tissue growth, since a state of saturated growth is unlikely to occur due to geometric and time constraints in a living organ. In addition, the distribution of the saturated residual stress increases the overall stress, and therefore acts contrary to the well-known observation that the residual stress reduces the stress in the wall under external loading. Perhaps it is for this reason, that the residual loop stress distribution at the saturation limit has not been observed in experiments on arteries or the heart. We also proved that our theory agrees with that of Goriely and Amar (2007) under the same assumptions of symmetric growth and deformation. However, for general elastic deformation, the updated growth tensor using our approach becomes a function of elastic deformation from all previous loadings. Given the wide range of growth and remodelling processes in the living tissues, our theory is an important step forward in studying the time-responses of live organs to external factors soft tissues.

Appendix: Thermodynamics compatibility
We now demonstrate that our constitutive approach for soft tissue growth satisfies the second law of thermodynamics. In other words, the system has a non-decreasing entropy 1 .
We consider the constitutive law (5) and (28) that includes the growth, but we drop the index k and the incompressibility constrain to be more general. For an elastic tissue with growth, the mass density 0 ( , t) in the reference configuration B 0 changes with time, and can be computed from the push back of the mass density ( , t) in the current configuration B t , i.e.
where J = det( t ) . Time derivative of 0 ( , t) gives where r is the mass source per unit volume, is the velocity vector and is the mass flux across the boundary B t .
Growth (or atrophy) is the result of mass production (or removal) and comes from a volumetric "source" r per unit volume in B t and the boundary influx m, together written where is the outward normal of B t . The linear momentum balance equation for the growing tissue then becomes where * = + (div ) , * = − ⊗ . Let k = 1 2 ⋅ be the kinetic energy per unit mass, the total kinetic energy of the system is and the rate of K is From the linear momentum equation (64) and the divergence theorem, we write the rate of force working P * as where = ∇ is the velocity gradient. The energy balance of the system including thermal effects is now where U = ∫ B t udv is the total internal energy, with u being the internal energy of the body per unit mass, G is the rate of energy supply to feed the growth (or the additional entropy production term), and Q is the rate of heat production defined as in which h is heart source per unit mass and q is the heat flux per unit area. (67) is then We now take so that last two terms cancel (and this not the only option, but is the simplest). We could think of G as consisting of bulk growth energy supply g (per unit mass) and energy flux (per unit area): thus, so locally we would have Epstein and Maugin (2000) included in their energy balance a flux term, which is essentially (k + u) , that we have not included here, but again this does not influence the construction of constitutive laws unless one brings in second-gradient effects (which considerably complex the theory). The same applies to the various dissipative terms they include in the energy balance and entropy production.

3
The local form of the energy balance (69) can now be expressed as Next define the free energy per unit mass by where T > 0 is the absolute temperature, and is the entropy per unit mass. Then, since the free energy is a function of density 0 , temperature T and deformation tensor, we can write (70) as With this we associate the entropy equation-in global form this is where H (a capital version of ) is the entropy production per unit mass, and in local form, or Combining this with (71), and assuming that the elastic material is isothermal (T is constant), we arrive at where we have made use of (5) and (28), following Ogden (1984), to get The second law of thermodynamics requires that H ≥ 0 , this is known as the entropy production inequality (the entropy of an isolated system is non-decreasing-it evolves towards thermodynamics equilibrium). From (74), we can see that this is satisfied, as in general for a growing tissue we have Thus, our constitutive modelling of the growing tissue is thermodynamically compatible as (70) u = tr( * T ) + h − div .
On the other hand, if the living body is atrophic, i.e.
the entropy production inequality (76) still holds.