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 bio-environment. 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 re-structure 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 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, 2009); 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 long-standing 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 so-called 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 aswhere is the pure growth tensor, and 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 , then a commonly used growth tensor can be expressed as (Kerckhoffs et al. 2012)where , , 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 (Kerckhoffs et al. 2012; Driessen et al. 2004), or combined stress- and strain- driven 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, and , respectively. The idea of growth from the reference configuration is illustrated in Fig. 1, where a pure (compatible) growth takes an elastic body from to an (imaginary) grown and stress-free configuration as , assuming that the map exists and is differentiable. Then, after loading the body deforms into the current configuration with elastic deformation gradient .

Fig. 1

Compatible growth from the reference configuration, . A pure and compatible growth takes the elastic body from into , and is then deformed via the loading-induced elastic deformation into the current configuration . The total deformation gradient from to is denoted as

The overall deformation gradient isFollowing (Ogden 2003), for incompressible materials that undergo volumetric growth, we must have the connectionwhere is the free energy per unit volume in , is the elastic strain energy function per unit volume in , , and . Since due to the hypothesis [H1], we simply havewhere we have used , and .

The Cauchy stress in can be determined either through the total deformation from , or through the elastic deformation from ,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 is made compatible via an elastic deformation into with residually stress, before the loading-induced maps into . As shown in Fig. 2 and [H2], this process can also be viewed as a growth-induced residual deformation from to followed by .

Fig. 2

A pure growth from takes the body to an incompatible and stress-free state in configuration . The deformation assembles the body into a compatible, but residually stressed stated in configuration . Loading induced deformation then maps it onto the current configuration . The process can be equally viewed as from to followed by from to , or from to directly via the elastic deformation

We now focus on incompressible materials and choose as the new (evolving) reference configuration, with a strain energy function . The total elastic deformation gradient from to is

As before, the Cauchy stress in can be determined either from the reference configuration , where is the elastic strain energy function related to the configuration .

or from ,where the strain energy function from is a function of both and residual stress , and p, are the Lagrange multipliers to ensure incompressibility of the material. The identity () 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 , we obtain (9).

When evaluated in with and , 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 , 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 . As shown in Fig. 3, is reached after the body is deformed from the reference configuration .

Fig. 3

Growth from a stressed configuration following an elastic deformation . Then, the body undergoes a small incremental deformation from to . The process can be equally viewed as from to followed by from to , or from to directly via the elastic deformation

Let indicate this elastic deformation. Since the line elements of material transform according to , we have the incremental connection,If the stressed configuration is now chosen as the reference configuration then the right-hand side of (14) becomes , where is the value of in this configuration. Here, is the unit tensor before the initial configuration : . Since and hence is independent of the reference configuration, we obtain the connection (Ogden 1997)Now consider an incremental deformation caused by the pure growth and , which take the body into the incompatible configuration first and then assemble it into a compatible configuration . If the growth displacement is “small” for each in so that terms of order are negligible in comparison with those of order then we refer to as an incremental (growth) deformation from to .

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 , , , (15) can also be written asAs in Sect. 2.1.2, the total elastic deformation gradient can be equally achieved either along the top path ––, or the lower path ––. The physical interpretation of the lower path –– is that a pure growth , which can be derived from , maps the reference configuration into the grown and incompatible (stress free) configuration . is made compatible by and becomes with the residual stress . From the body is deformed to by the elastic deformation . The process of the top path can be viewed as from to followed by from to . So again, we havefrom which we obtainWe 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 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 and .

Since pure growth is independent of the elastic stretch (H1), the only difference between and is due to the rotation of , i.e.where . Since can be easily computed, given , is known.

The rest of the analysis simply follows the approach of Sect. 2.1.2. The Cauchy stress in isThe strain energy function with respect to satisfiesThe residual stress is simply the Cauchy stress evaluated in when the external loading is removed, i.e. , . Hence,Again, and both obey the corresponding equilibrium equations, and the appropriate boundary conditions. From these, we solve for and , and and .

Subsequent growth from stressed and grown configuration

Subsequent growth can occur for k steps of growth and deformation in , as shown in Fig. 4. In general, external loading may also change after every step. Now with an new incremental growth from , the total deformation gradient from to becomeswhere ,is the deformation gradient from to , andHere, is the total cumulative pure growth from to , computed from the previous cumulative pure growth , i.e.Again, to follow the path ––, we must first find the connection between and , so that (26) can be evaluated. The deformation gradient between and is . Let , then we have the connectionwhere is the rotation of .

Fig. 4

General roadmap after k steps of continuous growth and deformation. The deformation gradient from to is , the incremental growth from the loaded configuration is . The total cumulative pure growth tensor to is . is the grown, stress free but incompatible configuration, and is the compatible and residually stressed configuration

The procedure is now the same as before. So we give the general expressions. The Cauchy stress in iswhereSimilar to (22), the residual stress isAgain, and both obey the corresponding equilibrium equations and the appropriate boundary conditions. From these, we can solve for and .

At this point, all we need is an incremental growth law for .

Incremental growth law

If we consider continuous growth and deformation, then we can write , similarly, . For simplicity, we assume that the growth tensor is diagonal for a fibre-reinforced 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 be the incremental growth tensor, we writewhere () is the fibre structure of the material at time t, which is related to the fibre structure at viaand in (31) are the rates of change of the growth multipliers in the directions, respectively, which can be written as . Note , so represents the total cumulative growth vector of the elastic body with respect to the reference configuration .

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 writewhere and are the Cauchy stress and Eulerian strain tensors in , 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. = 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 choosewhere are the maximum values of the growth multipliers , and , are growth parameters. The choice of and in (33) will ensure the growth rate changes smoothly until the growth multiplier has reached its maximum value . The nonlinearity and the speed of the growth are governed by and , respectively.

The strain energy function

We choose the grown and stress-free configuration 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) in ,where , are invariants of , the right Cauchy–Green tensor with respect to , is the total elastic deformation from to ,

and .

and are angles between and , and and in the plane, respectively, and

the fibre direction in changes according toAnd a, b , are material parameters.

At there is no growth, overlaps with , so .

Numerical methods

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.

Fig. 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 (, , ) for local coordinates, where , , and 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 – plane

At the current time t, after k steps of growth.

The growth evolution equation (33) predicts the incremental growthwhere .

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 principlewhere , are the element ordered vectors of stress and strain tensors, u is the displacement vector, is the ordered matrix of tangent moduli tensor (see (Taylor 2014) for details of the definition of an ordered matrix), t is the traction, is the symmetric part of , is the traction on the force boundary , is the displacement vector on the displacement boundary . The computational algorithm for the volumetric growth is given below.

Test examples

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

Modelling G&R in myocardium

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 (, , ), 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 . Note thatLet the myofibre architecture be described by a “fibre–sheet–normal” system (, , ) in the reference configuration,  (Wang et al. 2013). Here, we assume that the fibre direction always lies in the – plane, the sheet direction is transmural, and the sheet–normal .

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 kPa, . We also impose a constant internal pressure of 12 mmHg from throughout the growth.

LV case 1: athletic heart and isotropic growth

Since this is isotropic growth, we have such that (31) becomesWe note in this special case, the fibre structure () in the updated stress-free configuration is the same as the initial fibre structure ().

For Growth-, 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 lawwhereHere is the Mandel stress tensor, , , is the second Piola–Kirchhoff stress tensor, and is a critical stress value. The parameters are chosen asFor Growth-, we have the growth tensor (40) and its evolution law defined in , and (33) and (34), withand similarly choosewith all other parameters being the same as in Growth-.

Figure 6 shows that both models display ventricular dilation and wall thickening.

Fig. 6

Contours of the growth factor in the Case 1 after 0, 1, 3, and 5 weeks growth. Above: Growth-. Blow: Growth-. Path 1 at the upper middle section of the LV is defined to compare the transmural variation of results, see text for details

However, the growth patterns are clearly different. In Growth-, most of the growth is focused on the apex area, which is opposite to that of the Growth- 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, and , along the mean fibre direction, calculated from the residually stressed configuration as shown in Fig. 4. The transmural distributions of and 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-, is positive from the endocardial surface, and decreases transmurally; but for Growth-, is negative from the endocardial surface, and increases transmurally. The trend of Growth- 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.

Fig. 7

Transmural distribution of the Eulerian residual strain (a), and stress (b), in the mean fibre direction along path 1 for Growth- (dash-dot) and Growth- (solid), at weeks 1, 3, and 5. The corresponding experimental measurement of from a canine heart mid-anterior wall by Costa et al. (1997) is also plotted in (a) for comparison

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 assumingFor Growth-, the growth tensor is defined in asThe strain-driven involution law is Göktepe et al. (2010a),where is the same as in (41), is the stretch along the direction of , i.e. , and obeyswhere is the criteria value of elastic stretch. The parameters are chosen to befollowing Göktepe et al. (2010a).

For Growth-, we have the incremental growth tensor defined in ,We use the same strain-driven law as (43), except the function is function of ,where is the criteria value of elastic stretch along . In this case, the fibre structure in either or 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 replacing .

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-, and more growth away from the apex area for Growth-. 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.

Fig. 9

Contours of the growth factor in the Case 2 after 0, 1, 3, and 5 weeks growth. Above: Growth-. Below: Growth-. In both cases, the LV cavity size increases more obviously than the wall thickness

The residual strain and stress 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 conclusion are seen in (Göktepe et al. 2010b; Lee et al. 2015; Klepach et al. 2012).

Fig. 10

As in Fig. 7, but for the Example 2

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 bywhere L is the LV length, and , denote the inner and outer radii of each layer (), respectively.

In terms of cylindrical polar coordinates (r, , z), the geometry of the deformed configuration is given bywhere is the axial stretch. Similar to the incremental growth used by Goriely et al. (2008), we choose the incremental growth tensor in to beand the growth tensor for to bewith a simplified stress-driven growth involution law as in (Genet et al. 2015),andwhere and 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 of is the same as the growth tensor in .

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 aswhere is chosen to be 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 and stress for the 4-layer cylindrical model at different growth times, using Growth- and Growth- 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).

Fig. 11

Transmural distributions of a residual hoop strain , and b residual hoop stress , at growth time t =0.3, 1, 3, 5, for the 4-layer cylindrical model using Growth- (solid) and Growth- (dash-dot) approaches. Both models reach the same limit as time goes to infinite (), when an uniform growth is reached everywhere

Fig. 12

Transmural distributions of a the growth factor for Growth- (solid) and Growth- (dash-dot), (b) time history of the maximum residual hoop stress for Growth- (red solid with circle) and Growth- (black dash-dot with triangle), at growth time t=0.3, 1, 3, 5, for the 4-layer cylindrical model

Note that in both models, increases monotonically with time. However, the trend of stress/strain is different; it increases monotonically with time for Growth-, whereas for Growth- it reaches a local maximum first before decreasing. When growth is saturated, Growth- and Growth- 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 layersAfter the growth, the jth layer is deformed with the middle perimeter and the radius , where indicates the inner radius of the jth layer. Similarly, the middle radius of the next layer is where indicates the inner radius of the th layer. Material continuity (compatibility) requires thatSince , (59) can be rearranged to readIf we consider infinitely many layers, i.e. taking the limit of 0, then (60) is simplywhere 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 asIf , the growth will induce a smaller residual hoop stress in the jth layer compared with the th layer. If across the whole wall, then the hoop stress distribution will be from negative (compression) to positive (stretching). Likewise, if , the growth will induce a greater residual hoop stress in the jth layer compared with the th layer. If 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 across the whole wall initially in the 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 approach, we found 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 , giving sufficient time, all the material points eventually reach the same saturated growth in both approaches, so we havewhich leads to 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-, and the previous approach, Growth-, 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- model is supported by experimental observations. The other differences we noticed are that in the Growth- 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- approach, where the residual stress or strain increase monotonically with time until the set growth limit is reached, as shown in Fig. 8.

Fig. 8

The time history of the maximum residual hoop stress , at growth time t=1-5 weeks in case 1 for Growth- (red solid with circle) and Growth- (black dash-dot with triangle)

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 ( 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 , which is the same solution of the constant growth studied by Goriely et al. (2008) for .

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 from essentially plays the same role as the conventionally used growth tensor in (2). In other words, once we know , all the classical theories of G&R follow. However, it is impossible to define an evolution law for prior to deformation, since elements of 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 . 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 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 used by Goriely and Amar (2007)(i.e. “” in their Fig.2), is the same as in our approach. Note the rotation of is an identity tensor due to its diagonality (i.e. “ )) and “” is essentially the mapping from to in Fig. 4, using our notation this is . The connection of this mapping to is simply . Since all the tensors are diagonal and orders can be swapped, we can show thatIn 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 cannot be simplified to 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.