A mixed trigger volumetric growth law for cylindrical deformation in stressed configurations

Abstract

Soft tissue growth is crucial across various physiological applications, with mathematical modelling playing a pivotal role in understanding the underlying processes. The volumetric growth theory serves as a commonly used mathematical framework in this context. Our previous research on volumetric growth theory primarily concentrated on defining the incremental growth tensor in loaded and stressed configurations, revealing that this approach closely aligns with experimental observations of residual hoop stress distribution. However, given the assumptions employed, the approach has limitations in accurately predicting the growth timeline. In this work, we address these issues by incorporating the effect of initial residual strain and introducing a new mixed trigger growth evolution law. In this growth law, we do not use growth saturation as an upper limit, as this assumption cannot represent many physiological conditions. Instead, we propose that growth in soft tissues leads to a new equilibrium state. To illustrate this idea, we introduce a growth incompatibility function, denoted as $\mathcal{I}$ . We establish the analytical relationship between $\mathcal{I}$ and the opening angle in a simplified cylindrical geometry resembling the structure of the heart or arteries. We put forth a revised growth law that is both stress and incompatibility driven/Our results show that by using this mixed trigger growth law, tissues will not grow indefinitely. Instead, a stress-driven homeostasis incompatibility state will be reached. In addition, by accounting for the initial opening angle in the model, we can accurately trace the growth history of the heart, aligning with experimental data obtained from measuring the opening angle in young pigs from birth to maturity.

1. Introduction

Mathematical modelling of soft tissues plays a pivotal role in understanding growth and remodelling processes in physiological systems, such as the heart. This field has been actively researched for several decades, focusing on comprehending and simulating the dynamic growth and remodelling occurring within living tissues and organs, which hold significant relevance in both physiological and pathological contexts. Two prominent theoretical frameworks have emerged. The constrained mixture theory, pioneered by Humphrey and Rajagopal [1], views tissues as composed of multiple constituents with distinct natural configurations that can continuously change. Models rooted in this theory have been instrumental in investigating the growth and remodelling of arterial structures, aneurysms, and the heart [2]. However, accurately tracking the evolving natural configurations of all constituents remains challenging, especially in the absence of comprehensive experimental data.

An alternative framework, the volumetric growth theory, assumes a single reference configuration for all constituents. It decomposes the deformation gradient tensor into a growth tensor and an elastic tensor, effectively separating pure growth from elastic deformation. This concept, introduced by Rodriguez et al. [3], has been widely adopted in various applications, including Taber’s work on modelling enlarged arteries in 1998 [4] and the research conducted by Göktepe et al. [5] pertaining to the heart. Nevertheless, a key limitation lies in most models defining the growth tensor and its evolution law in the reference configuration, despite the fact that growth is an ongoing process within the current loaded configuration. Recent studies have endeavoured to model growth in the current configuration without releasing residual stresses but often relied on simplifications such as material symmetries [6]. A recent framework that accommodates arbitrary geometries and growth patterns fully evaluated in the loaded configuration was introduced by Zhuan and Luo [7]. This approach better captures the dynamic tissue adaptation to changing biomechanical environments and produces qualitative agreements with experimental observations in terms of the transmural residual stress distribution.

In our previous research, we identified that defining the growth tensor in stressed configurations allowed us to produce a qualitatively correct transmural distribution of residual strain within the heart and arteries. However, even with the implementation of the deformed growth tensor in the loaded configuration, traditional growth evolution laws primarily focusing on stress- or strain-driven mechanisms proved inadequate in capturing the time history of residual strain that is consistent with experimental observations. Our work shed light on the limitations of traditional stress-driven growth laws. These laws lead the system toward a state of saturated growth, as suggested by Moulton and Goriely [8]. However, this saturation limit is unlikely to be achieved, given the dynamic constraints imposed by the geometry and time frame of living organs.

In addition, most of the previous models ignored the initial residual strain induced by the heart formation. This omission led to an overestimation of blood pressure within arteries, resulting in discrepancies between estimated fibre growth trends and experimental data [9]. Furthermore, the embryonic heart tube initially carries substantial residual stresses, stemming not from growth but rather from the heart’s morphogenesis. Experimental studies by Taber et al. [10] and Taber [11] observed these phenomena during chick embryonic development, where heart membranes fused to form a single cardiac tube composed of three layers: the myocardium, the cardiac jelly, and the endocardium. The myocardium, composed of contractile cells, serves as the primitive heart’s only contractile element. These observations laid the foundation for exploring the intricate relationship between cardiac function, growth triggers, residual stress, and active stress.

The remodelling process in the heart is predominantly triggered by local stress conditions, primarily attributed to changes in cell shape rather than mere growth. The absence of normal cavity pressure leads to minimal growth during bending, as noted by Butler [12], and during in vivo looping, there is limited evidence to suggest that differential growth significantly contributes to the process, as indicated by Stalsberg and DeHaan [13].

A natural young heart undergoes a growth process as it matures, as defined by Sedehi and Ashley [14]. This enlargement of the heart is a physiological adaptation that occurs in response to increased pressure and thickening of its walls, ultimately reaching a state of growth equilibrium.

At the microscopic level, it has been documented that both residual strain and total stress/strain significantly influence the size of cardiomyocytes [15]. This influence is often characterized by a proportional increase in both the length and width of cardiomyocytes [16]. On a macroscopic scale, these cellular changes culminate in an increase in cardiac mass [17].

Pressure-induced changes drive ventricular wall thickening, while residual stress, stemming from growth incompatibilities, acts as a limiting factor for total stress. The interplay between residual stress and subsequent (Cauchy) stress reshapes the heart’s growth trajectory, guiding it toward a state of maturity and stability. This ideal state features both low residual stress and full heart functionality, ensuring the heart’s optimal performance.

Furthermore, it is essential to consider the presence of residual stress within the passive embryonic heart, as documented by Taber et al. [10]. The dorsal myocardium maintains longitudinal residual tensile stress through stretching or sustained contraction, leading to the bending of the cardiac tube with the dorsal myocardium located at the inner curvature. Even after the dorsal myocardium has resorbed, and the heart separates from the embryo, it retains its looped configuration, indicating a remodelling process to mitigate bending-induced stresses. This looping is a fundamental aspect of cardiac morphogenesis, involving biomechanical forces, even though the precise mechanisms remain incompletely understood.

From a growth mechanism perspective, looping introduces the initial opening angle, which triggers initial growth incompatibility and residual stress. This initial residual strain alters the growth path, influencing stress distribution and cell function during heart development. The implications of these experimental observations suggest the substantial presence of initial residual stress, originating from the heart’s formation. Subsequent growth mechanisms operate to reduce residual strains, guiding the heart toward an equilibrium state [18, 19].

Building on the concept of defining growth laws in loaded configurations, as introduced by Zhuan and Luo [7] in our study on volumetric growth, this paper introduces a novel mixed trigger growth law that incorporates both Cauchy stress and growth incompatibility as catalysts for growth. To elucidate this concept, we initially focus on a specific model system–the cylindrical geometry, relevant to the left ventricles and blood vessels.

To substantiate the qualitative predictions of our theory, we compare our findings with longitudinal measurements of the opening angles in the short-axis cut of the left ventricle in porcine hearts, tracked from birth. The congruence between the experimental trends in the opening angle over time and our theoretical predictions underscores the validity of our approach.

2. Kinematics of finite growth

Let $\mathbf{X}$ and $\mathbf{x}$ be the position vectors of a material point in the reference and current configurations, ${\mathcal{B}}_0$ and ${\mathcal{B}}_t$ , respectively. The pure stress-driven concept of volumetric growth from the reference configuration is illustrated in the lower path of Figure 1, where the total deformation gradient transforms a body from ${\mathcal{B}}_0$ to a grown and stressed configuration ${\mathcal{B}}_t$ . Assuming that the mapping $\mathbf{x}\left({\mathbf{X}}_0\right)$ exists and is differentiable, this overall deformation gradient is

$$\mathbf{A}={\mathbf{F}}_E{\mathbf{F}}_G={\mathbf{F}}_e{\mathbf{F}}_{\tau }{\mathbf{F}}_g,$$

(1)

where ${\mathbf{F}}_G$ and ${\mathbf{F}}_g$ are the total and incremental pure growth tensors, respectively, from ${\mathcal{B}}_0$ , ${\mathbf{F}}_{\tau }$ is the residual elastic deformation due to growth, ${\mathbf{F}}_e$ is the elastic deformation due to loading, and ${\mathbf{F}}_E={\mathbf{F}}_e{\mathbf{F}}_{\tau }$ is the total elastic deformation. Although in the first step as in Figure 1, we have ${\mathbf{F}}_G={\mathbf{F}}_g$ , in the subsequent growth, these two become different.

Figure 1.

Growth from a stressed configuration $\overline{\mathcal{B}}$ following a general deformation gradient ${\mathbf{A}}_1={\overline{\mathbf{F}}}_e$ . Then the body undergoes a small incremental deformation ${\overline{\mathbf{F}}}_g$ , which is made compatible by ${\overline{\mathbf{F}}}_{\tau }$ from ${\overline{\mathcal{B}}}_g$ to ${\mathcal{B}}_t$ . The result can be equally achieved following the lower path as ${\mathbf{F}}_g$ from ${\mathcal{B}}_0$ to ${\mathcal{B}}_g$ followed by ${\mathbf{F}}_{\tau }$ from ${\mathcal{B}}_g$ to ${\mathcal{B}}_{\tau }$ , and then ${\mathbf{F}}_e$ from ${\mathcal{B}}_{\tau }$ to ${\mathcal{B}}_t$ .

It is noted that in the pure stress-driven approach one often uses the form $\mathbf{A}={\mathbf{F}}_E{\mathbf{F}}_G$ only, because it is either assumed that growth does not induce residual strain, i.e. ${\mathbf{F}}_{\tau }=\mathbf{I}$ , as in Menzel and Kuhl [20], or that ${\mathbf{F}}_{\tau }$ is incorporated into ${\mathbf{F}}_E$ as here and in Taber and Eggers [21] and Goriely [6]. In this paper, since we explicitly address the residual deformation due to growth, we assume that normally ${\mathbf{F}}_{\tau }\ne \mathbf{I}$ , and write out equation (1)₂. This allows the pure growth-induced deformation ${\mathbf{F}}_{\tau }$ to be explicitly modelled.

We now briefly describe the theory that enables growth law to be defined in a loaded configuration [7]. As shown by the upper path in Figure 1, we denote such a configuration as $\overline{\mathcal{B}}$ . $\overline{\mathcal{B}}$ is reached following a deformation ${\mathbf{A}}_1$ from the stress-free configuration ${\mathcal{B}}_0$ .

Consider an incremental growth ${\overline{\mathbf{F}}}_g$ from the loaded configuration $\overline{\mathcal{B}}$ to ${\overline{\mathcal{B}}}_g$ , which is made compatible by ${\overline{\mathbf{F}}}_{\tau }$ (Note ${\overline{\mathbf{F}}}_{\tau }$ is not the same as the residual strain of the body that occurs when loading is removed), before arriving at ${\mathcal{B}}_t$ . Since both the upper and lower paths start from ${\mathcal{B}}_0$ and end at ${\mathcal{B}}_t$ , the results of these two paths in Figure 1 are equivalent, so we can also express the total deformation gradient as

$$\mathbf{A}={\overline{\mathbf{F}}}_{\tau }{\overline{\mathbf{F}}}_g{\overline{\mathbf{F}}}_e,$$

(2)

However, the significance of following the upper path is so that we can define an incremental growth tensor ${\overline{\mathbf{F}}}_g$ from the loaded configuration $\overline{\mathcal{B}}$ as in Zhuan and Luo [7]. If one follows the lower path in equation (1), then the incremental growth tensor ${\mathbf{F}}_g$ is always defined in a stress-free configuration.

Since we assume that pure growth is independent of the elastic stretch, the only difference between ${\overline{\mathbf{F}}}_g$ and ${\mathbf{F}}_g$ is due to the rotation of ${\mathbf{A}}_1$ , i.e.,

$${\mathbf{F}}_g={\mathbf{R}}_1^{\mathrm{T}}{\overline{\mathbf{F}}}_g{\mathbf{R}}_1,$$

(3)

where ${\mathbf{R}}_1={\mathbf{V}}_1^{-1}{\mathbf{A}}_1$ and ${\mathbf{V}}_1={\mathbf{A}}_1{{\mathbf{R}}_1}^{\mathrm{T}}.$ By making use of equation (3), we are able to obtain ${\mathbf{F}}_g$ and then proceed along the traditional growth and remodelling computational path from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ via ${\mathcal{B}}_g$ . In other words, the Cauchy stress in ${\mathcal{B}}_t$ is

$$\mathit{\sigma }={\mathbf{F}}_E\frac{\partial W_g\left({\mathbf{F}}_E\right)}{\partial {\mathbf{F}}_E}-p\mathbf{I},$$

(4)

where ${\mathbf{F}}_E={\mathbf{F}}_e{\mathbf{F}}_{\tau }$ . The strain energy function $W_g$ with respect to ${\mathcal{B}}_g$ is computed from the strain energy function $W_0$ with respect to ${\mathcal{B}}_0$ as

$$W_g\left({\mathbf{F}}_e\right)=J_g^{-1}{W}_0(\mathbf{A}),$$

(5)

where we have used $J_g=det{\mathbf{F}}_g$ and $J_e=1$ .

The residual stress $\mathit{\tau }$ is simply the Cauchy stress evaluated in ${\mathcal{B}}_{\tau }$ when the external loading is removed, i.e., ${\mathbf{F}}_e=\mathbf{I}$ , ${\mathbf{F}}_E={\mathbf{F}}_{\tau }$ ,

$$\mathit{\tau }={\mathbf{F}}_{\tau }\frac{\partial W_g\left({\mathbf{F}}_{\tau }\right)}{\partial {\mathbf{F}}_{\tau }}-p\mathbf{I}.$$

(6)

In other words, both $\mathit{\tau }$ and $\mathit{\sigma }$ obey the corresponding equilibrium equations and the appropriate boundary conditions. From these, we solve for ${\mathbf{F}}_{\tau }$ and ${\mathbf{F}}_e$ , and $\mathit{\tau }$ and $\mathit{\sigma }$ .

Subsequent growth can be computed after many steps of incremental growth and deformation, as shown in Figure 2, so $\overline{\mathcal{B}}$ becomes the ${\overline{\mathcal{B}}}_1$ , ${\overline{\mathcal{B}}}_2$ … and ${\overline{\mathcal{B}}}_n$ , ${\overline{\mathcal{B}}}_g$ becomes ${\overline{\mathcal{B}}}_{g_1}$ , ${\overline{\mathcal{B}}}_{g_2}$ , … and ${\overline{\mathcal{B}}}_{g_n}$ , ${\mathcal{B}}_G$ becomes ${\mathcal{B}}_{G_1}$ , ${\mathcal{B}}_{G_2}$ , … and ${\mathcal{B}}_G$ . Notice the external loading ${\overline{\mathbf{F}}}_{e_i}$ may change between each step. After $n\mathrm{th}$ steps of growth, A can be written as

$$\mathbf{A}={\overline{\mathbf{F}}}_{{\tau }_n}{\overline{\mathbf{F}}}_{g_n}{\mathbf{A}}_n={\mathbf{F}}_E{\mathbf{F}}_G,$$

(7)

where

$${\mathbf{A}}_n={\overline{\mathbf{F}}}_{e_n}...{\overline{\mathbf{F}}}_{{\tau }_2}{\overline{\mathbf{F}}}_{g_2}{\overline{\mathbf{F}}}_{e_2}{\overline{\mathbf{F}}}_{{\tau }_1}{\overline{\mathbf{F}}}_{g_1}{\mathbf{A}}_1,$$

and ${\mathbf{F}}_G$ is the total growth tensor with respect to ${\mathcal{B}}_0$ :

$${\mathbf{F}}_G={\mathbf{F}}_{g_n}...{\mathbf{F}}_{g_2}{\mathbf{F}}_{g_1},$$

(8)

with each ${\mathbf{F}}_{g_i}$ , $i=1,2,\dots ,n$ rotated from ${\overline{\mathbf{F}}}_{g_i}$ following equation (3), except ${\mathbf{A}}_1$ in equation (3) is replaced by the total deformation gradient ${\mathbf{A}}_i$ prior to the $i\mathrm{th}$ growth [7]. ${\mathbf{F}}_E$ is the total elastic deformation with respect to the grown configuration ${\mathcal{B}}_G$ :

$${\mathbf{F}}_E={\mathbf{F}}_e{\mathbf{F}}_{\tau }={\mathbf{AF}}_G^{-1}.$$

(9)

Figure 2.

General roadmap after $n$ steps of continuous growth and deformation, the deformation gradient from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ is $\mathbf{A}$ , the deformation gradient from the ${ℬ}_0$ to the loaded configuration ${\overline{\mathcal{B}}}_n$ is ${\mathbf{A}}_n$ . From ${\overline{\mathcal{B}}}_n$ we define the incremental growth ${\overline{\mathbf{F}}}_{g_n}$ , and the deformation is made compatible by ${\overline{\mathbf{F}}}_{{\tau }_n}$ . The total cumulative pure growth tensor ${\mathcal{B}}_0$ to ${ℬ}_G$ is ${\mathbf{F}}_G$ . ${\mathcal{B}}_G$ is the grown, stress-free but incompatible configuration, and ${\mathcal{B}}_{\tau }$ is the compatible and residually stressed configuration.

For a more general computational algorithm, please refer to Zhuan and Luo [7]. (Although ${\mathbf{F}}_{e_k}$ and ${\mathbf{F}}_{{\tau }_k}$ in equations (23) and (24) of Zhuan and Luo [7] should be ${\mathbf{F}}_{E_k}$ and ${\mathbf{F}}_{G_k}$ .).

3. Growth induced incompatibility in a cylindrical geometry

3.1. Concept of incompatibility index

To derive an appropriate growth law from the loaded configuration, we now restrict our analysis to a cylindrical model featuring isotropic deformation. Then ${\mathcal{B}}_0$ is described as

$$R_i\le R\le R_o,0\le \mathrm{\Theta }\le 2\pi ,0\le Z\le L.$$

(10)

Since the elastic energy of a growing isotropic material only depends on the symmetric part of the growth tensor [6, 22]. Hence, we only need to consider the case when ${\mathbf{F}}_g$ is symmetric. This means that ${\mathbf{F}}_g$ is no longer a two-field tensor, it is a one-field tensor with both legs defined in ${\mathcal{B}}_0$ (Figure 1). Likewise, ${\mathbf{F}}_G$ is fully also defined in ${\mathcal{B}}_0$ , and ${\overline{\mathbf{F}}}_g$ is fully defined in $\overline{\mathcal{B}}$ .

For simplicity, we further assume that the growth tensor is diagonal. This assumption is widely accepted by the community for the left ventricle and arteries (23). For a planar growth, the total growth tensor from equation (8) in matrix representation is

$$\left({\mathbf{F}}_G\right)=\mathrm{diag}({\vartheta }_R,{\vartheta }_{\mathrm{\Theta }},1),$$

(11)

where ${\vartheta }_R$ and ${\vartheta }_{\mathrm{\Theta }}$ are the cumulative growth multipliers in the $R$ and $\mathrm{\Theta }$ directions. Moreover, we assume that the cylindrical model geometry remains cylinder after deformation and growth, then both bases in ${\mathcal{B}}_0$ and $\overline{\mathcal{B}}$ are the same, and we have

$${\mathbf{F}}_g={\overline{\mathbf{F}}}_g,$$

(12)

This deformation also suggests that the growth parameters only spatially depend on the radial location, i.e.,

$${\vartheta }_R={\vartheta }_R(R),{\vartheta }_{\mathrm{\Theta }}={\vartheta }_{\mathrm{\Theta }}(R).$$

After deformation, we write $\mathbf{x}=\mathbf{\chi }(\mathbf{X},t)$ in ${\mathcal{B}}_t$

$$r_i\le r\le r_o,0\le \theta \le 2\pi ,0\le z\le l,$$

(13)

as

$$r=r(R),\theta =\mathrm{\Theta },z=\lambda Z.$$

(14)

where $\lambda$ is a constant pre-stretch parameter. Since the increase of the volume increase due to the growth is $J_G=det{\mathbf{F}}_G$ , we can write

$$r=\sqrt{\frac{{\int }_{R_i}^{R}2J_{G}R\mathrm{d}R}{\lambda }+r_i^2},$$

(15)

The total deformation is given as

$$\mathbf{A}=r^{\prime }{\mathbf{e}}_r\otimes {\mathbf{E}}_R+\frac{r}{R}{\mathbf{e}}_{\theta }\otimes {\mathbf{E}}_{\mathrm{\Theta }}+\lambda {\mathbf{e}}_z\otimes {\mathbf{E}}_Z,$$

(16)

where $\{{\mathbf{E}}_R,{\mathbf{E}}_{\mathrm{\Theta }},{\mathbf{E}}_Z\}$ and $\{{\mathbf{e}}_r,{\mathbf{e}}_{\theta },{\mathbf{e}}_z\}$ are the coordinate basis in ${\mathcal{B}}_0$ and ${\mathcal{B}}_t$ , respectively.

The material incompressibility requires that $det\left({\mathbf{F}}_E\right)=1$ , so we have

$$det\mathbf{A}=det\left({\mathbf{F}}_E{\mathbf{F}}_G\right)={\vartheta }_R{\vartheta }_{\mathrm{\Theta }}=J_G.$$

(17)

From equations (16) and (17), we have

$$\frac{r^{\prime }r}{R}\lambda ={\vartheta }_R{\vartheta }_{\mathrm{\Theta }}.$$

(18)

Substituting equation (18) into equation (15), we write

$$r^2=r_i^2+2{\int }_{R_i}^R\frac{1}{\lambda }{\vartheta }_R{\vartheta }_{\mathrm{\Theta }}R\mathrm{d}R.$$

(19)

To evaluate the residual strain ${\mathbf{F}}_{\tau }$ in ${\mathcal{B}}_{\tau }$ , we remove loading and let ${\mathbf{F}}_e=\mathbf{I}$ , then ${\mathbf{F}}_{\tau }={\mathbf{AF}}_G^{-1}$ , i.e.,

$$\left({\mathbf{F}}_{\tau }\right)=\mathrm{diag}(\frac{r^{\prime }}{{\vartheta }_R},\frac{r}{R{\vartheta }_{\mathrm{\Theta }}},\lambda ).$$

(20)

Now if the body is totally compatible then

$$F_{\tau }=I.$$

(21)

In other words,

$$\{\begin{array}{c}\frac{r^{\prime }}{{\vartheta }_R}=1,\hfill \\ \multicolumn{1}{c}{\frac{r}{R{\vartheta }_{\mathrm{\Theta }}}=1,}\\ \multicolumn{1}{c}{\lambda =1.}\end{array}$$

(22)

Rewriting equation (18) as

$$\frac{r^{\prime }}{{\vartheta }_R}=\frac{R}{r}{\vartheta }_{\mathrm{\Theta }},$$

(23)

combining with equation (22)₂, we see that equation (22)₁ automatically holds. From equations (22)₂ and (19), we obtain

$${\left({\vartheta }_{\mathrm{\Theta }}R\right)}^2=r_i^2+2{\int }_{R_i}^{R_o}{\vartheta }_R{\vartheta }_{\mathrm{\Theta }}R\mathrm{d}R.$$

(24)

Differentiate equation (24) with respect to $R$ , we have

$${\vartheta }_{\mathrm{\Theta }}^{\prime }=\frac{{\vartheta }_R-{\vartheta }_{\mathrm{\Theta }}}{R}.$$

(25)

Equation (25) can be viewed as the compatibility condition. Therefore, we introduce a growth incompatibility index as

$$\mathcal{I}=\frac{{\vartheta }_R-{\vartheta }_{\mathrm{\Theta }}}{R}-{\vartheta }_{\mathrm{\Theta }}^{\prime }.$$

(26)

We can see that $\mathcal{I}$ measures the total residual strain at ${\mathcal{B}}_{\tau }$ . $\mathcal{I}=0$ when the growth induced deformation is compatible and ${\mathbf{F}}_{\tau }=\mathbf{I}.$ If $\mathcal{I}>0$ , the local residual strain (stress) is positive, if $\mathcal{I}<0$ , the local residual strain (stress) is negative. Notice that $\mathcal{I}$ is a function of $R$ , i.e. it is a distribution.

3.2. A mixed-trigger growth law in a stressed configuration

It has traditionally been posited that the law governing growth evolution is driven by either stress or strain, in conjunction with a limiting function that prevents unrestricted growth. By defining incremental growth from a stressed configuration as ${\overline{\mathbf{F}}}_g=\mathrm{d}\mathit{\vartheta }$ , where $\mathit{\vartheta }=({\vartheta }_r,{\vartheta }_{\theta },{\vartheta }_z)$ , known as growth multipliers, Zhuan and Luo [7] investigated the growth of a left ventricle. Their modelling revealed that the distribution of transmural residual stress aligns with experimental observations when applying the current framework. Conversely, utilizing a pure stress-driven approach yields results that diverge from these observations.

However, the growth law used by Zhuan and Luo [7]

$$\stackrel{·}{\mathit{\vartheta }}=l(\mathit{\vartheta })\varphi (\mathit{\sigma }),$$

(27)

still involves a similar limiting function as in Göktepe et al. [5], i.e.,

$$l({\vartheta }_j)=\frac{1}{{\alpha }_l}{\left[\frac{{\vartheta }_j^{\max }-{\vartheta }_j}{{\vartheta }_j^{\max }-1}\right]}^{{\beta }_l},j=R,\mathrm{\Theta },Z,$$

where ${\vartheta }_j^{\max }$ are the maximum values of the growth multipliers, and the stress-driven trigger is defined as

$$\mathit{\varphi }(\mathit{\sigma })=\{\begin{array}{cc}{[\mathrm{tr}(\mathit{\sigma })-{\sigma }^{\mathrm{crit}}]}^{{\beta }_{\varphi }}\hfill & \mathrm{tr}(\mathit{\sigma })>{\sigma }^{\mathrm{crit}},\hfill \\ \multicolumn{1}{c}{0}& \mathrm{tr}(\mathit{\sigma })\le {\sigma }^{\mathrm{crit}}.\hfill \end{array}$$

(28)

where ${\sigma }^{\mathrm{crit}}$ is the critical value of Cauchy stress. ${\alpha }_l$ , ${\beta }_l$ , and ${\beta }_{\varphi }$ are genetic related growth parameters.

The introduction of the limiting function is primarily for mathematical convenience, as opposed to being grounded in experimental observations. Thus, in this work, we propose a new mixed trigger growth law. Our assumption here is that growth occurs in the current stressed configuration and is intricately driven by the interactions among stress, residual strain (or the incompatibility index), as well as the genetic growth. To this end, we modify the growth evolution law (27) to be

$$\stackrel{·}{\mathit{\vartheta }}=\psi (\mathcal{I})\varphi (\mathit{\sigma }),$$

(29)

where a non-negative function

$$\psi (\mathcal{I})=\frac{1}{{\alpha }_I}|1-\mathcal{I}/{\mathcal{I}}^{\mathrm{h}}{|}^{{\beta }_I},$$

(30)

is used to replace the limiting function. ${\mathcal{I}}^{\mathrm{h}}$ is a constant, representing the maximum local value of $\mathcal{I}$ at homeostasis, while ${\alpha }_I$ and ${\beta }_I$ are growth parameters. It is evident that when $\mathcal{I}=0$ , i.e., ${\mathbf{F}}_{\tau }=\mathbf{I}$ , growth is entirely driven by stress. However, what inhibits growth beyond the critical stress level? Inspired by experimental observations of growing hearts [18], we propose that the deforming body exhibits a tolerance threshold known as homeostatic incompatibility, which is often nonzero for mature organs. This mechanism effectively balances stress distribution when loading is applied either from within or outside cylindrical structures, such as the heart or blood vessels.

4. Thermodynamics considerations

To ensure that the new growth law is rational from the thermodynamics perspective, we invoke the Clausius-Duhem inequality. For a homogeneous temperature field from ${\mathcal{B}}_0$ to ${\mathcal{B}}_t$ , this inequality can be written as [6]

$$\mathbf{M}:\mathcal{G}+\overline{h}\ge 0,$$

(31)

where

$$\mathcal{G}={\stackrel{·}{\mathbf{F}}}_G{\mathbf{F}}_G^{-1}=\stackrel{·}{\mathbf{\Theta }}$$

(32)

is the rate of growth defined in the stress-free configuration ${\mathcal{B}}_0$ , $\mathbf{\Theta }=({\theta }_R,{\theta }_{\mathrm{\Theta }},{\theta }_Z)$ , $\mathbf{M}={\mathbf{F}}_E^{\mathrm{T}}\mathit{\sigma }{\mathbf{F}}_E^{-\mathrm{T}}$ is the Mandel stress defined in the virtual configuration ${\mathcal{B}}_G$ (Figure 2), and $-\overline{h}$ is the non-compliant entropy contribution to the process. Note, we have used $J_E=1$ in the expressions. In our cylindrical model, since $\mathbf{M}$ and $\mathit{\sigma }$ are coaxial, $\mathbf{M}=\mathit{\sigma }$ . This implies that the growth rate $\mathcal{G}$ should be a function of $\mathbf{M}$ or $\mathit{\sigma }$ to satisfy the inequality (31). In our cylindrical model, this function is chosen according to equation (29), given that $\mathcal{G}=\stackrel{·}{\mathbf{\Theta }}=\stackrel{·}{\vartheta },$ as shown in equation (12). Therefore, our proposed growth law aligns with thermodynamic principles.

Indeed, while thermodynamic considerations offer essential guidelines, they cannot prescribe a specific form for a growth law, particularly due to our limited understanding of $\overline{h}$ [5]. Since incremental growth invariably occurs within a stressed configuration, it is physically intuitive to suggest that such growth is driven by the actual (Cauchy) stress experienced by the tissues.

5. Application for porcine heart growth

We now apply our theory to the growth of porcine hearts. For simplicity, we consider a cylindrical heart model characterised by transversely isotropic growth. This model accounts for the proportional increase in cardiomyocyte length and width, leading to thickening of the left ventricular wall and an expansion of the chamber size, as observed by Hunter and Chien [16] and Pluim et al. [17].

Accordingly, the incremental growth tensor can be represented by a single growth multiplier $\vartheta$ as

$${\overline{\mathbf{F}}}_g/\mathrm{dt}=\stackrel{·}{\mathit{\vartheta }}=\stackrel{·}{\vartheta }\mathbf{I}+(1-\stackrel{·}{\vartheta }){\mathbf{e}}_z\otimes {\mathbf{e}}_z.$$

(33)

where

$${\stackrel{·}{\vartheta }}_r={\stackrel{·}{\vartheta }}_{\theta }=\stackrel{·}{\vartheta },{\vartheta }_z=1.$$

5.1. Opening angle measurements of growing porcine hearts

The open angle provides a simple quantitative measure of the residual strain distribution for incorporation into analytic and continuum mechanics models. The ventricular opening angle is defined as the angle between the two radial lines connecting the centre of the chamber, versus the centrelines of the cut wall (Figure 3). Here, porcine cardiac tissue was harvested from animals aged 1 day, 14 days, 7 months and 3 years. A heart from each age is presented for comparison in Figure 4(a)–(d). The youngest heart was harvested from an animal who died due to hypoxia either during, or immediately after, unassisted farrowing (birth). All other tissue was harvested from animals slaughtered for the food industry. All hearts appeared healthy and fully developed, with consistent tissue colour and no obvious abnormalities. A scalpel (blade size 11; Swann Morton, UK) was used to cut the hearts in the short axis, through the equatorial region. A single, 6 - 8 mm thick slice, was dissected from each heart and just submerged in cardioplegic fluid, which counteracted the gravitational forces and meant the slice was in a no-load state. The slice was then photographed from above, before cutting the left ventricle free wall radially, opposite the interventricular septum. A second photograph was then taken less than 30 s after the cut (Figure 4(d)). Figure 3 presents two photos of a slice, one before and one after the cut, demonstrating measurement of the open angle.

Figure 3.

Definition of opening angle. Left: Drawing of a slice in the no-load state. The chamber centre of this configuration is estimated by constructing a diameter through the epicardial edge of the cut and the most distant epicardial point. Right: Same slice in the opening angle state. The opening angle is defined as the angle between the two radial lines connecting the centre of the ventricular chamber and the centre lines of the walls at the cut edges.

Figure 4.

Harvested porcine hearts of four age groups, namely 1 day old (a), 14 days old (b), 7 months old (c), 3 years old (d). Scale bar = 20 mm.

For each group of three porcine hearts, the averaged measurements of the opening angles at the four time instants are $93.2^o\pm 8.6^o$ , $75.1^o\pm 12.6^o$ , $33.4^o\pm 6.2^o$ and $30.7^o\pm 4.1^o$ , respectively.

5.2. Computation of the incompatibility index based on opening angle measurement

For time steps $1,2,\dots ,n$ , each incremental growth tensor of the heart model is computed from equation (33). Furthermore, $\mathbf{A}$ , ${\mathbf{F}}_E$ , and ${\mathbf{F}}_G$ are computed from equations (7)–(9), respectively.

Material incompressibility requires that $det{\mathbf{F}}_E=1$ . The invariant $I_1$ of right Cauchy-Green deformation tensor ${\mathbf{F}}_E^{\mathrm{T}}{\mathbf{F}}_E$ is

$$I_1={\lambda }_{E_r}^2+{\lambda }_{E_{\theta }}^2+{\lambda }^2,$$

where

$${\lambda }_{E_r}=\frac{r^{\prime }}{\vartheta },{\lambda }_{E_{\theta }}=\frac{r}{R\vartheta }.$$

The components of the Cauchy stress tensor in ${\mathcal{B}}_t$ are then

$$\begin{array}{c}{\sigma }_{\mathrm{rr}}=-p+\frac{\partial W_G}{\partial I_1}{\lambda }_{E_r}^2,\hfill \\ {\sigma }_{\theta \theta }=-p+\frac{\partial W_G}{\partial I_1}{\lambda }_{E_{\theta }}^2,\hfill \\ {\sigma }_{\mathrm{zz}}=-p+\frac{\partial W_G}{\partial I_1}{\lambda }^2,\hfill \end{array}$$

where $W_G$ is the isotropic part of the invariant-based constitutive law for the myocardium in ${\mathcal{B}}_G$ developed by Holzapfel and Ogden [24].

In the absence of body forces, the stress components ${\sigma }_{\mathrm{rr}}$ and ${\sigma }_{\theta \theta }$ in ${\mathcal{B}}_t$ must satisfy the equilibrium equation $\mathrm{div}\mathit{\sigma }=0$ , i.e.,

$$\frac{d{\sigma }_{\mathrm{rr}}}{\mathrm{dr}}+\frac{{\sigma }_{\mathrm{rr}}-{\sigma }_{\theta \theta }}{r}=0.$$

Given a heart pressure $P$ and assuming zero axial force, the governing equations become

$${\int }_{r_i}^{r_o}\frac{{\sigma }_{\mathrm{rr}}-{\sigma }_{\theta \theta }}{r}=P,$$

(34)

$$\pi {\int }_{r_i}^{r_o}(2{\sigma }_{\mathrm{zz}}-{\sigma }_{\mathrm{rr}}-{\sigma }_{\theta \theta })r=0.$$

(35)

By solving equations (29) and (35), we determine $r_i$ and $z$ in ${\mathcal{B}}_t$ and, consequently, the Cauchy stress. Then, at each growth step $(1,2,\dots ,n)$ , applyiing the growth law (29) and solving equations (34) and (35) again with $P=0$ , we can calculate the residual stress tensor $\mathit{\tau }$ in ${\mathcal{B}}_{\tau }$ .

To estimate the opening angle caused by growth at time $t$ , we introduce a radial cut at ${\mathcal{B}}_{\tau }$ denoted by coordinates $(r,\theta ,z)$ , and let the cylindrical model relax into a grown and incompatible configuration with coordinates $(\rho ,\varphi ,\xi )$ and an opening angle $\alpha$ , such that

$${\rho }_i\le \rho \le {\rho }_o,\frac{\alpha }{2}\le \varphi \le 2\pi -\frac{\alpha }{2},0\le \xi \le L.$$

It is important to note that, in general, this configuration does not exactly correspond to ${\mathcal{B}}_G$ . However, if we assume that a single radial cut can release all the residual strain–a simplification as observed by Zhuan and Luo [25]–then it can be considered a mechanical proxy for ${\mathcal{B}}_G$ . Henceforth, the configuration following an opening angle cut will be referred to as ${\mathcal{B}}_G$ .

The total deformation tensor $\mathbf{A}$ from ${\mathcal{B}}_0$ to ${\mathcal{B}}_{\tau }$ is

$$\mathbf{A}={\rho }^{\prime }{\mathbf{e}}_{\rho }\otimes {\mathbf{E}}_R+\frac{\rho k}{R}{\mathbf{e}}_{\varphi }\otimes {\mathbf{E}}_{\mathrm{\Theta }}+{\mathbf{e}}_{\xi }\otimes {\mathbf{E}}_Z$$

where

$$\rho =\sqrt{k{\int }_{R_i}^{R}2J_{G}R\mathrm{d}R+{\rho }_i^2},$$

(36)

$J_G=\det (\mathbf{A})$ , and

$$k=2\pi /(2\pi -\alpha ).$$

(37)

Since the loading is removed, we set ${\mathbf{F}}_e=1$ . Consequently, the residual elastic deformation $F_{\tau }$ from ${\mathcal{B}}_G$ to ${\mathcal{B}}_{\tau }$ is computed using equation (9). Therefore, the corresponding residual stress components in ${\mathcal{B}}_{\tau }$ must satisfy the equilibrium equation $\mathrm{div}\mathit{\tau }=0$ . By solving equations (34) and (35) under the boundary condition $P=0$ in ${\mathcal{B}}_{\tau }$ , we determine ${\rho }_i$ . Subsequently, we calculate $k$ and the opening angle $\alpha$ from equations (36) and (37).

The corresponding incompatibility distribution is then derived as

$$\mathcal{I}=\frac{{\rho }^{\prime }}{R}(1-\frac{1}{k})=\frac{k}{\rho J_G}(1-\frac{1}{k})=\frac{k-1}{\rho J_G}.$$

(38)

Therefore, given a measured opening angle (or $k_h$ ) from a mature heart, the homeostatic incompatibility index is chosen as

$${\mathcal{I}}^h=\mathrm{sgn}\left(\frac{k_h-1}{{\rho }_m{J}_{G_h}}\right)\underset{\rho ={\rho }_m}{max}\left|\frac{k_h-1}{\rho J_{G_h}}\right|,$$

(39)

where $J_{G_h}$ is the determinant of the total deformation gradient at homeostasis. Assuming that $n$ steps of growth occurred at time $t$ , and $m$ steps of growth occurred when homeostatic incompatibility was reached, we have

$$J_G=J_{g_n}\dots J_{g_2}{J}_{g_1},\mathrm{and}{J}_{G_h}=J_{g_m}\dots J_{g_2}{J}_{g_1}.$$

(40)

5.3. Heart growth with an initial opening angle

We now consider a heart model with an initial opening angle ${\alpha }_0$ , equivalent to starting from the incompatible configuration ${\mathcal{B}}_G$ :

$${\rho }_i\le \rho \le {\rho }_o,\frac{{\alpha }_0}{2}\le \varphi \le 2\pi -\frac{{\alpha }_0}{2},0\le \xi \le L.$$

Using the cylindrical coordinates $(\rho ,\varphi ,\xi )$ for ${\mathcal{B}}_G$ , and $\{r,\theta ,z\}$ for the corresponding intact configuration ${\mathcal{B}}_{\tau }$ , we have

$$r=r(\rho ),\theta =k_0(\varphi -\frac{{\alpha }_0}{2}),z=\lambda \xi ,$$

where $k_0=2\pi /(2\pi -{\alpha }_0).$

By closing the initial opening angle, the incompressibility condition becomes

$$r=\sqrt{\frac{{\rho }^2-{\rho }_i^2}{k_0\lambda }+r_i^2}.$$

(41)

Solving the equilibrium equation (34) with $P=0$ , we obtain the inner radius $r_i$ of the closed cylinder in ${\mathcal{B}}_{\tau }$ . The residual strain tensor, due to the closure of the opening angle, is given by

$${\mathbf{F}}_{{\tau }_0}=r^{\prime }{\mathbf{e}}_r\otimes {\mathbf{E}}_{\rho }+\frac{r{k}_0}{\rho }{\mathbf{e}}_{\theta }\otimes {\mathbf{E}}_{\varphi }+\lambda {\mathbf{e}}_z\otimes {\mathbf{E}}_{\xi },$$

(42)

following [26].

We denote ${\mathbf{F}}_{G_0}$ as the underlying growth required to take an imaged stress-free and intact reference configuration ${\mathcal{B}}_0$ to ${\mathcal{B}}_G$ . Using the coordinate $(R,\mathrm{\Theta },Z)$ for ${\mathcal{B}}_0$ , we have

$$R_i\le R\le R_o,0\le \mathrm{\Theta }\le 2\pi ,0\le Z\le l.$$

Thus, for ${\mathcal{B}}_G$ , we can write

$$\rho =\rho (R),\varphi =\frac{1}{k_0}\mathrm{\Theta }+\frac{{\alpha }_0}{2},\xi =Z.$$

In other words, the deformation gradient from ${\mathcal{B}}_0$ to ${\mathcal{B}}_G$ is

$${\mathbf{F}}_{G_0}=\frac{\partial \rho }{\partial R}{\mathbf{e}}_{\rho }\otimes {\mathbf{E}}_R+\frac{\rho }{R}\frac{\partial \varphi }{\partial \mathrm{\Theta }}{\mathbf{e}}_{\varphi }\otimes {\mathbf{E}}_{\mathrm{\Theta }}+\frac{\partial \xi }{\partial Z}{\mathbf{e}}_{\xi }\otimes {\mathbf{E}}_Z,$$

(43)

From which we can deduce ${\mathbf{F}}_{G_0}=\mathrm{diag}({\vartheta }_R,{\vartheta }_{\mathrm{\Theta }},1)$ as

$$\begin{array}{c}\multicolumn{1}{c}{{\vartheta }_R=\frac{\partial \rho }{\partial R}={\rho }^{\prime },}\\ \multicolumn{1}{c}{{\vartheta }_{\mathrm{\Theta }}=\frac{\rho }{R}\frac{\partial \varphi }{\partial \mathrm{\Theta }}=\frac{\rho }{R{k}_0},}\\ \multicolumn{1}{c}{{\vartheta }_Z=\frac{\partial \xi }{\partial Z}=1,}\end{array}$$

and

$$\rho =\sqrt{k_0\underset{R_i}{\overset{R}{\int }}2J_{G_0}R\mathrm{d}R+{\rho }_i^2},$$

(44)

where $J_{G_0}=\det \left({\mathbf{F}}_{G_0}\right)$ . Notably, as the deformation is from a closed configuration to one with an opening angle, equation (44) is different to equation (41).

Accordingly, the initial incompatibility distribution is

$${\mathcal{I}}^0=\frac{k_0-1}{\rho J_{G_0}}.$$

(45)

We now solve the problem following the same procedure as before, except that the total growth tensor at time $t$ becomes

$${\mathbf{F}}_G={\mathbf{F}}_{g_n}...{\mathbf{F}}_{g_2}{\mathbf{F}}_{g_1}{\mathbf{F}}_{G_0}.$$

(46)

It is important to note that although $J_{G_0}$ is generally unknown, its exact value is not required for our calculations. This is because, in equation (30), $J_{G_0}$ is effectively cancelled out in the ratio $\mathcal{I}/{\mathcal{I}}^h$ , as both terms $\mathcal{I}$ and ${\mathcal{I}}^h$ have the factor $J_{G_0}$ .

6. Results

In all the simulations that follow, the parameters are chosen as follows :

$${\alpha }_l=2,{\beta }_l={\beta }_{\varphi }=1,{\sigma }^{\mathrm{crit}}=0.0012\mathrm{MPa},P=140\mathrm{mmHg},$$

in line with our previous paper [7]. In addition, for equation (30), we choose ${\alpha }_I=2,{\beta }_I=1,$ $1.5$ , and $2$ . ${\mathcal{I}}^0$ is computed from the opening angle measured at time zero, with ${\mathcal{I}}^{\mathrm{h}}=0.78$ estimated from measurements at $t=3$ years.

6.1. Pure stress-driven growth law

We first apply the pure stress-driven growth law, as defined in equation (29), to the heart model. The transmural distribution of residual hoop stress, ${\tau }_{\theta \theta }$ , is depicted in Figure 5(a). This distribution qualitatively agrees with the experimental observations by Costa et al. [27]. In our previous paper [7], we demonstrated that this residual stress distribution results from defining the incremental growth tensor in the stressed configuration.

Figure 5.

Transmural distributions of the hoop residual stress ${\tau }_{\theta \theta }$ using the pure stress-driven growth law without opening angle (a) and with the opening angle (b). Blue-squared line represents the residual stress estimated by Wang et al. [28] based on opening-angle experiment on heart by Costa et al. [27]

However, the slope of the distribution begins to increase, reaches its peak at week 14, and then gradually declines to zero as growth saturates. In other words, all material points converge toward the same upper limit, in line with predictions by Goriely [6]. At this upper limit, the growth incompatibility index drops to zero across the myocardidum. This mechanism accounts for the eventual reduction of residual stress to zero when applying the pure stress-driven growth law, a result that diverges from experimental findings.

The results from employing the same model, albeit with an initial opening angle of ${\alpha }_0=\frac{1}{2}\pi$ , are presented in Figure 5(b). The distributions again show similarities. The maximum hoop residual stress value increases from 5.12 kPa at week 0 to 7.30 kPa by week 7, peaking at 11.22 kPa around week 14, before it gradually diminishes, returning to its original value as time approaches infinity.

The overall trend is depicted more clearly in Figure 6, showcasing the time history of the opening angle, both with and without the initial opening angle considered. In the case lacking an initial opening angle, the opening angle increases steadily from the start until it reaches a maximum. It then diminishes as all material points meet their growth upper limits. With the initial opening angle factored in, the calculated opening angle initially surges to a peak of $\frac{1}{2}\pi$ , before gradually reverting to its initial state over time.

Figure 6.

Time history of the opening angle (rad) using the pure stress-driven growth law without opening angle (a) and with the opening angle (b)

In summary, introducing an initial opening angle alters the time-dependent behaviour of residual stress, making it similar to scenarios without an initial angle but changing the long-term upper limit. This adjustment may mirror real-world situations where disease-induced growth prompts an organ to revert to its previous homeostatic state once health is restored, as suggested by Humphrey and Schwartz [29]. However, such a growth pattern does not account for the healthy development of a young porcine heart, which is the focus of our study. Ramasubramanian et al. [30] observed that in a chick embryo heart, the opening angle decreases from its initial value, eventually stabilising at a mature equilibrium with a smaller opening angle due to heart wall thickening, which reduces initial residual stress.

6.2. Mixed trigger growth law

We now apply the mixed trigger growth law, (30), for the cylindrical heart model, with and without the initial opening angle. The transmural distributions of the hoop residual stress are shown in Figure 7 for different time instants. As shown in Figure 7(a), the maximum residual stress increases from 0 to a peak value of 1.03 kPa at week 9. Afterward, it gradually decreases, approaching nearly zero by the tenth week, even without imposing an upper limit on the growth multiplier as in equation (27). This behaviour results from the mixed trigger growth law, wherein both growth incompatibility and Cauchy stress affect circumferential growth.

Figure 7.

Transmural distributions of the hoop residual stress ${\tau }_{\theta \theta }$ at weeks 1–35 using the mixed trigger growth law without (a) and with the opening angle (b). The squared data are obtained by Guan et al. [2].

The transmural distribution of Cauchy stress, following the application of inner pressure, diminishes from the endocardium to the epicardium, a trend that is inverse to that of residual stress, aligning with our simulations and those of others [2]. This disparity promotes more rapid growth from the inner surface and slower growth from the outer surface. Conversely, the growth incompatibility, indicative of residual strain, shifts from negative at the endocardium to positive at the epicardium. As a result, incompressibility prompts more accelerated growth from the outer surface and decelerated growth from the inner surface. Therefore, Cauchy stress and incompatibility function as competing triggers, culminating in a reduction of residual strain (or stress) over time.

The transmural distributions of residual stress, influenced by the initial opening angle, are illustrated in Figure 7(b). Here, the maximum residual stress diminishes steadily from 6.03 kPa at week 1 to 0.51 kPa by week 21, ultimately reaching a state of equilibrium. This pattern is attributed to the interplay between growth triggers: growth incompatibility and Cauchy stress. Initially, the significant growth incompatibility induced by the initial opening angle seeks resolution, leading to a decline in residual stress. As this process unfolds, the effect of Cauchy stress remains relatively consistent. Over time, an equilibrium is established between Cauchy stress and growth incompatibility, resulting in a stable distribution of residual stress.

The temporal evolution of opening angles under the mixed trigger growth law is depicted in Figure 8(b), eventually stabilising at a significant yet smaller angle. This behaviour closely aligns with our experimental measurements in porcine hearts, where the opening angles measured at various stages showed a gradual decrease, mirroring the computed trend that also reaches a stable limit as indicated by experimental observations. The consistent decrease in the opening angle is due to the thickening of the heart wall, aimed at reducing the initial residual stress. As the heart matures, the opening angle stabilises at a lower level, sufficient for maintaining heart function. This reduced equilibrium residual strain state has also been observed by other researchers [6, 31].

Figure 8.

The time history of the opening angle (rad) using the mixed trigger growth law, with no initial opening angle (a) and with initial opening angle (b). The measured opening angles of young pigs’ hearts from 1 day to 30 weeks and 3 years are shown as symbols and error bars in (b). It is noted that homeostasis is reached in all cases after 30 weeks. Furthermore, it appears that a value of ${\beta }_I=1.5$ provides the best description of the experiments. Although values of 1 and 2 also yield qualitatively good results.

7. Discussion

The core objective of this study is to introduce a new mixed trigger growth law that captures incremental growth from the loaded configuration, aiming for a qualitative alignment with experimental observations. Our investigation is currently limited to cylindrical models and isotropic materials, and we have applied this model to the study of porcine heart growth. For the cylindrical deformation, we have shown that the proposed new growth law satisfies thermodynamic conditions. We acknowledge that real heart geometry significantly differs from a cylindrical model, and that myocardial tissue is far from being an isotropic material. In this context, the remodelling process due to changes in fiber structures is not considered here.

In our previous work [2, 7], guided by stress-driven growth laws, we applied the general framework to more complex models, incorporating anisotropic growth in both cylindrical and more realistic three-dimensional heart models. This included modelling directional growth along fibres, sheets, and normals. However, this paper does not explore the complexities associated with the effects of fibres (myocyte and collagen fibres). To develop a new growth law and facilitate comparison with experimental measurements of heart growth, we opted for a simplified cylindrical model. This choice is in line with a well-established precedent, as the cylindrical heart model has been extensively utilised by numerous researchers, particularly in experimental studies [5, 10, 32–36]. Integrating fibre growth and remodelling into our model would necessitate more complex tensor rotations and the strain energy functions–areas we intend to explore in future studies.

The reason we assume that the z-direction is growth-free in our examples is due to the experimental measurements being limited to radial cuts in the middle wall of the heart. In the absence of specific data, the simplest assumption is that there is no growth in this direction. However, we acknowledge that this may not accurately reflect reality and have addressed this point in the Discussion section. The fact that our growth law has achieved qualitative agreement with experimental observations for the first time suggests that any potential growth in the z-direction does not play a dominant role.

We understand that growth processes are typically multiphysics in nature, as they involve the coupling of various phenomena such as mass transport, nutrient diffusion, and mechanical deformation. Consequently, the stimuli for growth can be diverse. In what is known as homeostatic growth, which occurs during tissue maintenance, stress acts as the direct stimulus for growth [37]. In contrast, in many cases, growth is propelled by nutrient diffusion. In our growth law, we introduce a tuning factors, ${\alpha }_I$ and ${\beta }_I$ , to partially regulate the growth rate. These parameters could be adapted in a nutrient-diffusion growth model, e.g., function (3) in Soleimani et al. [38], where nutrient concentration $C$ is utilised to regulate growth. Accordingly, we can integrate the stimuli of stress, residual strain, and nutrient diffusion to control the growth of living organs in future.

The validation of our model is underscored by its alignment with experimental data. Considering the theoretical nature of our work, the room for parameter variation is limited; the homeostatic value of the incompatibility index, ${\mathcal{I}}^h$ , was deduced from measurements, whilst other parameters were obtained from existing literature. To examine the sensitivity of our results to the value of ${\beta }_I$ in (30), which should depend on nutrients and other growth factors, we conducted simulations for ${\beta }_I=1,1.5,$ and $2$ , finding the model’s predictions to be robust. Similarly, altering ${\alpha }_I$ has a smaller effect than ${\beta }_I$ (not shown), primarily influencing the timing of equilibrium establishment. Notice we could also use the distribution of the incompatibility index for ${\mathcal{I}}^h$ , not a local maximum value. However, our simulations show that the effect is similar to changing the value of $\beta$ . Given that the configuration following the opening angle cut is only a proxy, we believe that choosing the maximum local value of the incompatibility index is simpler and equally rational.

Finally, we would like to discuss the differences between our methodology and those methodologies employing homeostatic surfaces as described in Lamm et al. [22] and Holthusen et al. [39]. In the homeostatic surface approach, a growth potential is introduced as a function of the conjugated driving forces and a set of material parameters, representing a hemostatic surface located within the principal stress space. The growth tensor is defined so as to minimize the amount of energy needed to reach this homeostatic state. In other words, the direction of growth is the gradient of the potential to the Mendal stress in this framework. In our approach, the growth is controlled by reducing either the deformation incompatibility or the trace of Cauchy stress. The key difference is that they defined the incremental growth tensor in the reference configuration, and therefore can consider the elastic energy and growth energy as decoupled. In our approach, the incremental growth is defined in the stressed configuration, so in general it is driven by the Cauchy stress as well as by previous growth and loading history [7].

Our current growth law is specifically tailored for cylindrical models, which have many clinical applications. However, there is a need to develop a more universal growth law for general deformation that can be applied across a variety of situations. An intriguing avenue for future research may be exploring the integration of homeostatic surface-driven growth with incremental growth laws defined within stressed configurations.

8. Conclusion

In this study, we adopted a new mixed-trigger approach using both stress and growth incompatibility to jointly drive the growth process. In addition, our work introduced the effect of the initial opening angle, a factor often overlooked in previous studies. We demonstrated that considering this initial condition is essential in the growth history of the newborn heart. To illustrate this concept, we applied our methodology to a simplified multi-layer cylindrical model, one that allows for inhomogeneous growth in residually stressed and loaded configurations. Subsequently, we compared our findings concerning residual stress distribution and opening angle with experimental measurements taken from healthy young porcine hearts. Furthermore, our results revealed that without considering the initial opening angle, the residual strain increases first and then reduces to a homeostatic limit as time extends to infinity, a trend unsupported by our experimental data. This is becuase the residual strain in the new formed heart is often at its maximum level. However, our mixed trigger growth law combined with the consideration of initial opening angle, yielded results that significantly aligned with experimental observations of growing porcine hearts. In essence, our estimated opening angle demonstrated a reduction from its initial value and decrease to an equilibrium over time. This trend not only agrees with our experimental measurements but also concurred with a well-documented phenomenon: the reduction of stress in the tissue wall under external loading as the newborn heart matures.