Critical roles of time-scales in soft tissue growth and remodeling.

Most soft biological tissues exhibit a remarkable ability to adapt to sustained changes in mechanical loads. These macroscale adaptations, resulting from mechanobiological cellular responses, are important determinants of physiological behaviors and thus clinical outcomes. Given the complexity of such adaptations, computational models can significantly increase our understanding of how contributions of different cell types or matrix constituents, and their rates of turnover and evolving properties, ultimately change the geometry and biomechanical behavior at the tissue level. In this paper, we examine relative roles of the rates of tissue responses and external loading and present a new rate-independent approach for modeling the evolution of soft tissue growth and remodeling. For illustrative purposes, we also present numerical results for arterial adaptations. In particular, we show that, for problems defined by particular characteristic times, this approximate theory captures well the predictions of a fully general constrained mixture theory at a fraction of the computational cost.

INTRODUCTION

As aptly stated in 1995 by Fung, “Every specialty in biomechanics begins with the study of constitutive relations,” that is, descriptors of material behavior for particular conditions of interest. Among his many contributions, Fung showed that although soft tissues tend to exhibit nonlinearly viscoelastic behaviors under certain conditions, they can often be regarded as pseudoelastic under physiological conditions. Indeed, it is for this reason that preconditioning is a fundamental part of most experimental protocols for studying the biomechanical behavior of soft tissues and the vast majority of constitutive relations for stress are based on hyperelasticity, not viscoelasticity. When considering the remarkable ability of soft tissues to respond to changing mechanical conditions, that is to grow and remodel, Fung stated further that “the scope of constitutive equations broadens: it should now include a mass-and-structure growth-stress relationship as well as a stress-strain relationship.” For the past 20+ years, a key specialty area of investigation in soft tissue mechanics has focused on developing and testing constitutive relations for growth and remodeling (G&R).

Of the different approaches available, we have advocated and employed a constrained mixture theory for soft tissue G&R. Briefly, this approach requires identification of three classes of constitutive relations: a hyperelastic descriptor of the mechanical behavior of each of the structurally significant constituents, as well as descriptors for mass density production rates and related survival functions for these same constituents. The mechanical properties and rates of production and removal of the individual constituents can depend, in part, on the state of stress/strain at which they were incorporated within the extant tissue, as well as on the current state of stress/strain. Hence, this approach is consistent with Fung's call for “growth-stress” relations. A special case of tissue maintenance exists when rates of production and removal balance perfectly while constituent's turnover in an unchanging mechanical state. Finally, the term “constrained” implies that all motions of each constituent correspond with those of the mixture even though each constituent can possess an individual natural (i.e., stress-free) configuration. This constrained mixture approach has proven useful in describing a host of evolving vascular conditions, including the development and resolution of cerebral vasospasms, mechano-adaptation of arteries to altered blood pressure and flow, arterial aging, the enlargement of aneurysms, the development of tissue engineered vascular grafts, and the maladaptation of vein grafts, as well as other cellular and tissue-level processes.

Because every cell and structurally significant constituent has a finite half-life and presumed memory associating its loss to the state of stress at which time the constituent was incorporated within the extant matrix, the classical constrained mixture theory uses a hereditary integral formulation similar to that of nonlinear viscoelasticity. Albeit motivated directly by the mechanobiology, this integral formulation can be expensive computationally, including the need to store all past states over which constituents were deposited. For this reason, there has been a search for suitable simplifications that preserve advantages of the mixture approach (e.g., the ability to account for material properties and rates of turnover inherent to the different constituents that constitute the tissue) while improving computational efficiency. One approach has been to introduce a temporal homogenization while another has been to derive an associated steady-state form that reveals the final (evolved) state, not unlike equivalence derivations in viscoelasticity that relate integral and rate forms or those that focus on long-term responses in relaxation and creep. In contrast, in this paper we consider time scales inherent to the rates of mechanical loading and G&R responses to determine conditions under which a rate-independent (“pseudoelastic”) theory can hold throughout G&R. In this sense, our current formulation is similar to concepts introduced by Fung to model tissues that exhibit viscoelastic behaviors using concepts of hyperelasticity. For purposes of illustration and application, we use this theory to simulate arterial responses to altered pressure, flow, and axial stretch.

RESULTS

Long-term, steady-state, tissue maintenance solution

The goal of this first example is to confirm that the rate-independent (pseudoelastic) G&R model derived in Sec. IV B can compute exactly the long-term response (Sec. IV C) of a thin-walled bilayered artery (Sec. IV A 1) when subjected to multiple external loads that are sustained for long periods. Consider, therefore, two different combinations of loading consisting of 1.15-fold increases in the applied pressure P and flowrate Q, each sustained following initial transients. To delineate better the responses to different stimuli, the loads are applied sequentially in the orders P(1) → Q(2) (first case) and Q(1) → P(2) (second case), each taking 21 days overall to reach steady values and time-shifted from the other by 14 days. See Fig. 1(a).

FIG. 1.

Predictions by the full constrained mixture model (first case, dashed; second case, dash-dotted) for the evolution of (c) medial and adventitial collagen (referential) mass density, (d) inner radius, (e) medial thickness, and (f) adventitial thickness, each normalized to original values, which result from two different cases of perturbations in loading [(a), solid lines] that cause (b) evolving stimulus functions (i.e., stresses different from homeostatic target values). Panels (g)–(i) show, separately, associated deviations in stress components from original homeostatic values. Note the two different time-delayed combinations of changes in pressure and flow (a). Shown, too, is the long-term, mechanobiologically equilibrated solution (solid square), which is the same for each of the two (same final) loading conditions. Note the perfect correspondence of the long-term steady-state solution computed with the present time-independent formulation and the full (hereditary) constrained mixture model.

Panel (b) shows the resulting/evolving stimulus function ϒ for (both medial M and adventitial A) collagen c as an example. Panels (c)–(f) in Fig. 1 show the case-specific evolving responses (for , a, h, and h, respectively) predicted by the full model of Sec. IV A for the different combinations of loads shifted over time. Finally, fold differences in shear, circumferential, and axial stresses from homeostatic (present in ) are shown separately in panels (g)–(i). Note, in particular, the complex responses that are obtained by initiating the different mechanical perturbations at different times. Indeed, such simulations illustrate the importance of modeling G&R because results are not always intuitive at first given the many parallel nonlinear processes. For example, an instantaneous increase in pressure would be expected first to increase luminal radius and decrease (isochorically) the wall thickness, whereas an infinitely slow increase in pressure might result in fully adapted (quasi-equilibrium) restorations of luminal radius and increases in wall thickness at each time. Actual G&R would be expected to fall within these extremes depending on rates and extents of loading and rates of matrix turnover.

Here, because the characteristic rate of G&R (k = 1/7 ≈ 0.143 days−1) is greater than the characteristic rate of change of the external loads (in this example, k ∼ 0.15/20 ≈ 0.075 days−1, recall Sec. IV C), both G&R processes start [Figs. 1(b)–1(f)] right after the first external load (either P or Q) increases [Fig. 1(a)]. Note, too, that the mechano-stimulus functions for collagen ϒ [Fig. 1(b)] and smooth muscle ϒ (not shown), which drive the G&R process, yield different short- and mid-term mass density [Fig. 1(c)] and geometric [Figs. 1(d)–1(f)] evolutions, yet regardless of the order in which the loads are applied, they yield a common long-term outcome at s = 105 days = 15s ≫ s = 7 days, when mechanobiological equilibrium is restored [mathematically described by , recall Eq. (31)] and the mass production rates balance perfectly the removal rates in the evolved and now unchanging configuration. Particularly interesting is the resetting of homeostatic stresses (τ ≠ τ, σ ≠ σ, and σ ≠ σ) at the new evolved state, which, yet satisfy the more general equilibrium condition .

Importantly, we can also directly compute this long-term, path-independent solution using the particularized formulation of Sec. IV B, which yields a single solution for the combination of external loads γ = P/P = ε = Q/Q = 1.15, here computed via a Newton–Raphson method with initial guess given by an ideal adaptation of the type , and , see Sec. IV E. This solution is shown by solid squares in Fig. 1, which reveals precise correspondence with the full, time-dependent constrained mixture model at long (fully adapted) times.

Instantaneously adapted, quasi-equilibrium, slow evolution

As discussed in Sec. IV C, rate-independent formulations derived in Sec. IV B are also valid for computing slow G&R if the characteristic rate k is much greater than the characteristic rate of change of the external stimuli k, satisfying then the quasi-steady-state relation k/k ≪ 1 at any G&R time s. For these particular situations, the arterial adaptation (occurring over a time scale s = 1/k) to different alterations in mechanical loading (occurring over a time scale s = 1/k ≫ 1/k = s) may be regarded as immediate at any G&R time s. We confirm this assessment in this example, quantifying at the same time how short the G&R characteristic time s should be with respect to the external loading characteristic time s so that the predictions given by both formulations, general and particularized, become (approximately) equal. We analyze three different cases, one for each type of external perturbation that stimulates the G&R response of the artery under study (i.e., altered inner pressure, flow rate, or axial stretch), showing the different adaptive processes that the artery undergoes for each.

Figure 2 shows three different G&R arterial responses (non-solid lines), associated with three different G&R characteristic times s = 1/k = 0.1 days, 1 day, or 10 days (with ), which are computed with the full model of Sec. II A for the particular increases in inner pressure γ(s) ≡ P(s)/P shown in Fig. 2(a). In order to assess the quasi-steady-state assumption for each case, these values of s are compared to the characteristic time of the external load application, which we estimate from Fig. 2(a) as the time taken for P/P → 1.15, namely s ∼ 10 days in this case. We also show in the same figure the single rate-independent solution of Sec. IV B (pseudoelastic, solid line) computed as a function of the time-dependent inner pressure ratio P(s)/P where time s simply plays the role of a simulation-driver parameter. In this case, we start the iterative solution procedure at each new time step using the converged solution at the previous time step.

FIG. 2.

Rate-independent (solid line) and rate-dependent (non-solid lines) evolutions computed, respectively, with the pseudoelastic and the full constrained mixture models, the latter with different characteristic G&R times s = {0.1, 1, 10} days, for an isolated increase in pressure with s ∼ 10 days. Shown are (a) prescribed load P/P from 1 to 1.15, (b) mechano-stimulus function ϒ, (c) referential mass densities of collagen , (d) relative inner radius a/a, (e) relative medial thickness h/h, and (f) relative adventitial thickness h/h. The final total wall thickness is h = 0.0494 mm (with 67% due to medial thickening and 33% due to adventitial thickening). Finally, shown too is the long-term mechanobiologically equilibrated solution (solid square), which reveals perfect correspondence of all three methods at the final adapted state. The scales are the same in Figs. 2, 3, and 4 to facilitate comparisons.

FIG. 3.

Similar to Fig. 2, except for an isolated increase in flow rate Q from 1 to 1.15. The final total wall thickness is h = 0.0407 mm (with 70% due to medial thickening and 30% due to adventitial thickening).

FIG. 4.

Similar to Fig. 2 except for an isolated increase in axial stretch λ from 1 to 1.1. The final total wall thickness is h = 0.0419 mm (with 68% due to medial thickening and 32% due to adventitial thickening).

Figure 2 shows that the prescribed increase in pressure provokes simultaneous changes in referential mass densities and [Fig. 2(c), with similar tendency for medial smooth muscle], inner radius a [Fig. 2(d)], and layer thicknesses h and h [Figs. 2(e) and 2(f)], driven by mechano-stimulus functions ϒ [Fig. 2(b) for collagen α = c, with similar tendency for smooth muscle], in four simulations (three rate-dependent based on different characteristic times s and one rate-independent). In particular, the evolution of inner radius, thicknesses, and masses for s/s ∼ 0.1/10 = 0.01, for which ϒ ≃ 1 [cf. Eq. (48)], is (in practice) indistinguishable from the rate-independent evolution, for which [cf. Eq. (31)]. For the other two simulations, s/s ∼ 1/10 = 0.1 and s/s ∼ 10/10 = 1, the mechano-stimulus function no longer satisfies the quasi-equilibrium condition ϒ ≃ 1 during early times. The evolution of luminal radius for these two cases initially separates from the mechanobiologically equilibrated solution, even though all remain close to the initial homeostatic value for our prescribed modest increase in pressure, i.e., a/a ≃ 1 = ε1∕3 where ε = Q/Q. Note, however, that the overall solution (including medial and adventitial thicknesses and constituent masses) given by the rate-independent formulation is still in good agreement with that of the full formulation. Finally, we see again that the rate-independent, all three rate-dependent, and the mechanobiologically equilibrated (solid square) simulations predict exactly the same long-term, fully equilibrated solution at s = 100 days, reaching a value (not shown) of the evolved total thickness h/h = 1.23 which is 7% greater than that for an ideal mechanoadaptation 1.15 = γε1∕3. This 7% difference between the final value of relative total thickness (1.23) and the ideal target (1.15) can be attributed to the relatively high content of elastin within the artery under study (33% overall), which is known to prevent a perfect adaptation since elastin does not turnover. Albeit not shown, the responses of shear, circumferential, and axial over-stresses allow one to understand the transient, short-term trends in Fig. 2, as, for example, in the case of comparable timescales (s/s ∼ 1). Indeed, the increase in pressure mainly provokes an initial increase in circumferential stress, hence the stimulus functions ϒ drive a growth process (with parallel increments of thickness and mass) until, eventually, stresses return close to normal values and .

Figure 3 shows the same type of analysis but for a particular increase in flow rate ε(s) ≡ Q(s)/Q [Fig. 3(a)]. We see that the rate-independent (pseudoelastic, solid line) formulation again provides a very good approximation to the quasi-static evolution predicted by the full model (case ) and exactly the same long-term, tissue-maintenance solution as the full model in all cases (, s = 100 days). If we compare the full model predictions for s/s ∼ 0.1 and s/s ∼ 1 to the single one given by the rate-independent model, we observe good agreements for the evolution of inner radius, but initially opposed tendencies for thicknesses and masses. We note that, at s = 100 days, and (not shown), indicative of a near but not fully mechanoadaptive solution again. In the case of comparable timescales (s/s ∼ 1), the increase in flow rate mainly provokes an increase in shear stress, and thus a decrease in the stimulus functions ϒ [cf. Eq. (11)], which initially attenuates matrix turnover consistent with nitric oxide slowing collagen production by smooth muscle cells (with parallel decrements of thickness and mass). This reduced thickness, along with the increase in luminal radius (also consistent with vasodilation with nitric oxide), provokes an increase in circumferential stress that, subsequently, drives a growth process via ϒ > 1. Finally, stresses closely return to normal values, such that .

Similar conclusions regarding quasi-equilibrium () and long-term (, s = 100 days) predictions are obtained for a single increase in the axial stretch ratio λ(s)/λ (in this case up to λ/λ = 1.10, Fig. 4). Regarding short- and mid-term evolutions, especially for s/s ∼ 1, we see that all geometric and mass variable responses separate from the quasi-equilibrium solution. This finding is consistent with prior results that mechanobiological adaptations are particularly sensitive to changes in axial length. Finally, the common mechanobiologically equilibrated state (not fully reached at s = 100 days for s/s ∼ 1) is such that a/a ≃ 1 = ε1∕3 and h/h = 1.05 > 1 = γε1∕3. In this case, for comparable timescales (s/s ∼ 1), the increase in axial stretch mainly provokes an increase in axial stress and thus stimulus functions ϒ, which initially drive a growth process (with parallel increments of thickness and mass). Yet, the substantial reduction in inner radius (because of the axial stretch) also provokes an increase in shear stress, which, subsequently, attenuates G&R via values ϒ < 1. Finally, because of the decreased thickness, circumferential stress increases slightly, driving G&R until the artery reaches a mechanobiologically equilibrated state associated with . Interestingly, the mechano-stimulus functions ϒ [panel (b) in Figs. 2–4 for the specific case of collagen] go through one extremum only (i.e., a maximum) in the case of an isolated increase in pressure, two extrema (i.e., a minimum and a maximum) in the case of an isolated increase in flow rate, and three extrema (i.e., a maximum, a minimum, and a last maximum) in the case of an isolated increase in axial stretch, which, in any case, equal the number of primary mechanical stimuli that are being stimulated sequentially.

Finally, Fig. 5 shows results similar to those in Fig. 2 except for three different degrees of increased luminal pressure γ(s) = P(s)/P for the slowest response examined in Fig. 2 (s/s ∼ 1). As it can be seen, increases in the perturbation in pressure from 1.15 to 1.30 to 1.45 fold results in progressively slower adaptations as expected. Nevertheless, in each case, the pseudoelastic model (solid curves) provides the same long-term predictions as the full model (dotted curves) and even provides reasonable predictions over the short-term in some metrics.

FIG. 5.

Rate-independent (solid lines) and rate-dependent (dotted lines) evolutions computed, respectively, with the pseudoelastic and the full constrained mixture models, the latter with a single characteristic G&R time s ∼ 10 days, for different increases in pressure with s ∼ 10 days. Shown are (a) prescribed loads P/P from 1 to 1.15, 1.30, and 1.45, (b) respective mechano-stimulus functions ϒ, (c) referential mass densities of collagen , (d) relative inner radius a/a, (e) relative medial thickness h/h, and (f) relative adventitial thickness h/h.

DISCUSSION

Biological growth (changes in mass) and remodeling (changes in microstructure) processes are necessarily time dependent due to the finite periods needed for the significant material to be synthesized, deposited, degraded, and/or reorganized. Moreover, rates of degradation appear to depend in part on the conditions present at the time of deposition (e.g., how the constituent was oriented or cross-linked), hence endowing tissues with a “material memory.” For this reason, hereditary integral formulations, as commonly used in viscoelasticity, have proven useful in describing and predicting evolving geometries, compositions, and properties of soft tissues under diverse conditions.

Nevertheless, such formulations can be computationally expensive and there is strong motivation to identify additional methods for analysis as well as the conditions for which such methods hold. In this paper, we presented a time-independent (“pseudoelastic”) formulation of arterial G&R to describe particular behaviors of otherwise time-dependent processes, which parallels Fung's use of pseudoelasticity to describe particular behaviors of otherwise viscoelastic tissues. We show that this approach, represented by a system of nonlinear algebraic equations for a bilayered, thin-walled artery, yields the same outcomes as a general constrained mixture model, represented by a system of nonlinear integro-differential equations, for both long-term steady state responses in which the external loads are sustained over time and quasi-equilibrium evolutions in which the external loads change slowly enough that the arterial wall can essentially adapt instantaneously to the given alterations at each G&R time s. In the latter case, the quasi-equilibrium formulation gives good qualitative results even for some cases in which rate-dependent effects remain, especially for isolated increases in blood pressure. This last observation, along with a marked reduction in computational time, makes the present pseudoelastic G&R formulation a good starting point for studying more complex situations that necessitate a general solution given by the full integral model. Importantly, the time needed to compute the rate-independent formulation (in evolution form) was ∼10 to 100 times less than that for the full model (when run in MATLAB on a 12 dual-core processor, and depending on specific values of s, which affect the time step of the integral formulation).

Obviously, even though (rate-independent) pseudoelasticity will never replace (rate-dependent) viscoelasticity, hyperelastic models describing the pseudoelastic behavior of viscoelastic tissues continue to play important roles during stages of experimental characterization, the computation of important responses, both analytically and numerically, the delineation of limiting responses, and in straightforward determinations of long-term outcomes following relaxation and creep. In parallel, time-independent G&R models as presented herein will never replace time-dependent G&R models, either fully integrated or kinematically motivated. We submit, however, that pseudoelastic G&R formulations can play parallel roles as in Fung's pseudoelasticity to better understand complex G&R of soft tissues, in general, and of arteries, in particular. Furthermore, because full and mechanobiologically equilibrated formulations predict the same long-term steady state outcomes, these simplified formulations can serve as efficient tools for constrained mixture modeling that includes studies of parameter sensitivity, uncertainty quantification, and optimization, which tend to be even more demanding computationally. Regarding material characterization of evolving tissues, the present formulation can be used, for example, to determine G&R-related material parameters by comparing original homeostatic and evolved homeostatic configurations [via, for example, Eqs. (39) and (40)], hence simplifying this important stage of constitutive modeling. Finally, this type of modeling G&R could also guide the process of engineering new tissues that react to artificially induced mechanical loading, as in Refs. 30–32.

For purposes of illustration, we included characteristic times (half-lives) for G&R that may be regarded unrealistic except in cases of extreme adaptations or pathologies (e.g., s ∼ 0.1 days), which we used as a more stringent test of the concept. Nevertheless, values of s should be assessed relative to a characteristic time of changes in the loads that stimulate G&R, namely s. Thus, a dimensionless number of the type s/s, comparing characteristic times of G&R and loading, is what actually determines the goodness of a quasi-static assumption. Regarding actual arterial behaviors, s may range from ∼10 to 100 days whereas changes in loading may range from seconds to many days, hence dimensionless numbers s/s of order unity, or lower, exist for s ∼10 to 100 days or longer. The lower the value of s/s, the better the pseudoelastic G&R approximation, as one would expect from viscoelasticity theory. In contrast, the greater the ratio s/s, the better the transient “elastic” (without G&R) approximation; see Example 2 in Ref. 13. Because s is material-dependent, it can change during certain conditions or diseases, hence each situation should be evaluated individually.

Timescales for G&R and mechanical loading can be very different in other situations, but analyses in terms of orders of magnitude are yet useful. Indeed, kinematic models of growth ultimately rely on the recognition of these two time scales, with evolution of the growth tensor F depending on the “growth (or) remodeling timescale.” Hence, growth laws may be expressed in kinematic models in terms of rate parameters that, alternatively, may depend on problem-specific growth metrics to prevent unlimited growth, with baseline values identified as rates of initial growth. Similarly, the so-called global growth approach also postulates evolution equations for the stress-free state of an artery, where, again, characteristic times associated with the growth process can be identified.

In conclusion, we emphasize that all of the illustrative results (Figs. 1–5) also depend on all of the particular constitutive assumptions, including functional forms (Secs. IV A and IV E) and prescribed values (Table I). Indeed, whereas we must continue to search for mechanobiologically appropriate and yet computationally tractable theoretical frameworks, the primary need—as noted by Fung decades ago—remains identification of the best constitutive relations, which for constrained mixture theories of G&R includes descriptors of mechanical behaviors, rates of mass production, mechanisms of incorporation within extant matrix, and rates of mass removal for each structurally significant constituent. Improved relations for complex pathologies representing maladaptation remain wanting. Theoretically motivated experiments are thus essential, which as we show herein should begin to focus more on the time scales over which loading and biological processes occur.

TABLE I.

Baseline material parameters for a mouse descending thoracic aorta. Best-fit values determined from Ref. 21, and adapted, for the specific examples performed in this work. In particular, values would typically be of order 1/70 days−1 under normal conditions, but we consider here a faster rate of turnover as a more stringent test for the pseudoelastic G&R framework. Superscripts e, m, and c denote elastin, smooth muscle, and collagen, with superscripts/subscripts r, θ, z, and d denoting radial, circumferential, axial, and symmetric diagonal directions. Subscripts M and A denote medial and adventitial layers whereas o denotes original homeostatic values. Subscripts σ and τ denote intramural and wall shear stresses, respectively. See Table II for definitions of other variables.

Artery mass density	ρ	1050 kg/m3
Medial mass fractions	[ϕMoe,ϕMom,θ,ϕMoc,z,ϕMoc,d]	[0.4714, 0.4714, 0.0381, 0.0190]
Adventitial mass fractions	[ϕAoe,ϕAoc,θ,ϕAoc,z,ϕAoc,d]	[0.0333, 0.0175, 0.3201, 0.6291]
Diagonal collagen orientation	α0	45.36°
Inner radius, thicknesses	[ao, hMo, hAo]	[0.6468, 0.0284, 0.0118] mm
Elastin material parameter	ce	114.5 kPa
Smooth muscle parameters	[c1m,c2m]	[401.0 kPa, 0.012]
Collagen parameters	[c1c,c2c]	[411.2 kPa, 5.5]
Deposition stretches	[Gre,Gθe,Gze]	[1/(1.9 × 1.6), 1.9, 1.6]
Deposition stretches	[Gθm=Gθc,Gzc,Gdc]	[1.071, 1.193, 1.192]
Maximum active stress	Tmax	100 kPa
Active stretch limits	[λM, λ0]	[1.1, 0.4]
Vasoactive parameters	[CB, CS]	0.8326× [1, 0.5]
Vasoactive rate parameter	kact	1/7 day−1
Mass production gains	[Kσm,Kτm,Kσc,Kτc]	[1.6, 2, 2, 2.5]
Mass removal rates	[kom,koc]	[1/7, 1/7] day−1

TABLE II.

Definition of volume-specific variables. MPR, Mass Production Rate; SEF, Strain Energy Function; α = {e, m, c}, Γ = {M, A}.

Mass of constituent α within layer Γ per unit current volume of layer Γ	ρΓα
Mass of constituent α within layer Γ per unit reference volume of layer Γ	ρΓRα
MPR of constituent α within layer Γ per unit reference volume of layer Γ	mΓRα
Nominal MPR of constituent α within layer Γ per unit reference volume of layer Γ	mΓNα
Volume-specific SEF of constituent α (defined at the constitutive level)	W^^α
SEF of constituent α within layer Γ per unit reference volume of layer Γ	WΓRα

METHODS

A bilayered constrained mixture model for arteries

First consider the kinematics and evolution equations of a constrained mixture model that governs G&R processes of an idealized bilayered artery that comprises a medial and an adventitial layer. Each layer is modelled as an independent constrained mixture having its own local variables that must satisfy kinematic compatibility constraints at the medial-adventitial interface. This time-dependent, integral-type formulation will be particularized in Sec. IV B to special cases for which the time dependence either vanishes or can be neglected, following ideas introduced in Ref. 13.

Geometry, mass fractions, and kinematics

Figure 6 illustrates an idealized bilayered cross-section of a normal elastic artery wherein the media (inner layer, since the normal intima is not expected to carry significant loads) is the functional layer and the adventitia (outer layer) bears increasingly more load as the pressure increases abruptly. Three illustrative in vivo configurations are shown at three respective G&R times (on orders of days to weeks). Let the (original) homeostatic configuration at time s = 0, namely κ = κ(0), serve as a reference. The configuration at the current G&R time s > 0 is denoted κ(s), while configurations at all intermediate times 0 ≤ τ ≤ s are denoted by κ(τ). Hereafter, for the sake of notational simplicity, we will omit the time dependence when it is not needed explicitly. We denote the inner (luminal) radius of the artery in a generic configuration as a, the medial thickness as h, and the adventitial thickness as h. The total wall thickness is thus h = h + h. The length of the artery is l, which is the same for both layers by compatibility of axial displacements at the medial-adventitial interface r = a + h, hence l = l ≡ l.

FIG. 6.

Evolving in vivo configurations of a bilayered arterial wall from time s = 0 (original homeostatic reference configuration κ(0) ≡ κ) to an arbitrary time s > 0 [current configuration κ(s)], showing the different deformation gradients F and F experienced by the medial and adventitial layers over time, in particular at times 0 ≤ τ ≤ s. Shown, also, are the different natural configurations for the different structurally significant constituents α (elastin “e,” smooth muscle “m,” and collagen “c”), which are deposited with separate but constant deposition stretches G within both layers at the indicated deposition times τ (except for the smooth muscle, present in media only), as well as geometric parameters a (luminal radius), h (medial thickness), and h (adventitial thickness) and luminal pressure P, which can also change over time.

We analyze the wall mechanics using a thin-walled approach in which variables within each layer are uniform, though different, consistent with effects of residual stress [Residual stresses tend to homogenize transmural distributions of stress (Refs. 16 and 17), thus rendering this assumption reasonable]. Hence, we separately compute the deformation gradient tensor for each layer, which quantifies motions between the original homeostatic reference configuration κ and a generic configuration κ. Using cylindrical coordinates X = {r, θ, z} to denote the position of a material point within either the media (M) or the adventitia (A), and assuming that cylindrical axes coincide with principal directions of strain, the respective layer-specific deformation gradients read and with stretches and where r = a + h enforces compatibility of radial (and circumferential) displacements at the medial-adventitial interface. Subscript o refers to the original homeostatic state.

The material composition is also layer-specific. In particular, we consider three primary types of load-bearing constituents (elastin-dominated “e,” smooth muscle “m,” and fibrillar collagen-dominated “c”) , with mass fractions satisfying in the media, and satisfying in the adventitia, with . Following a constrained mixture approach for each layer, the motion of each constituent (and cohort thereof) is constrained to equal the motion of the corresponding layer as a whole, as given by deformation gradients F and F in Eqs. (1) and (2). Each constituent and cohort, however, has its own evolving natural configuration , where τ denotes the instant at which that constituent mass was produced and deposited within the arterial wall (see Fig. 6). Without a loss of generality, we let the different constituents deposit in the media and adventitia with layer-independent deposition (symmetric) stretch tensors G, which means that the natural configuration of a constituent α at time τ, namely, , is common for both layers (see Fig. 6). In addition, we assume (in physiologic adaptations) constant and volume-preserving deposition stretch tensors G, with . Following the deformation path shown in Fig. 6, the deformation gradient experienced by constituent α deposited at time τ within layer Γ = M, A (when applicable) that survives to current G&R time s reads for smooth muscle and collagen, which are continuously produced and removed, as and for elastin, which is produced perinatally and thus at s ≤ 0 in our case of G&R in maturity, as since FΓ(τ = 0) = I.

Mass and strain energy evolutions

Smooth muscle cells and collagen fibers are produced continuously, hence their removal is usually accounted for in constrained mixtures models through mass survival (exponential decay) functions of the type (we use the form introduced in Ref. 13, but others are possible) where the rate “parameter” , with τ ≤ t ≤ s, is assumed (constitutively) to increase with respect to its original homeostatic constant value through with Δσ(t) quantifying any relative difference between a given scalar measure of intramural Cauchy stress (e.g., magnitude, principal invariant, maximum principal value) acting at time t at the tissue level, namely , and its corresponding homeostatic value , such that

Additionally, production of collagen and smooth muscle is governed constitutively by respective mass density production relations. Following Ref. 13 for the mass production rate of cohort α per unit reference volume of the mixture (in this case, within layer Γ), we have for both layers (see Tables I and II for nomenclature) where is an evolving nominal mass production rate including, importantly, the same function employed for mass removal and the referential mass density of constituent α within layer Γ (i.e., per unit reference volume of the respective layer Γ). This relation ensures balanced production and removal in homeostatic states, for which . Moreover, is a “stimulus function” that ultimately drives mass production to rates different from nominal depending on biochemomechanical stimuli; herein, we consider mechano-stimuli only. Over the years, we have found that a reasonable form for for arteries subjected to perturbed blood pressures and flows, linearized about the initial homeostatic state, reads where and are constituent- and layer-specific (constant) gain parameters, τ is the flow-induced shear stress acting over the endothelium, and Δτ = (τ – τ)/τ. In particular, at the original homeostatic state, where Δσ = 0 = Δτ.

Given constitutive equations for production and removal, the evolution of the layer-specific referential mass density of the cohort α is Assuming that the spatial mass density ρ of the arterial wall (i.e., its current mass per unit current volume of mixture) remains constant, the strain energy function of constituent α within layer Γ, also defined per unit reference volume of layer Γ, reads where is the volume-specific strain energy function of constituent α and is the right Cauchy–Green deformation tensor obtained from , which for smooth muscle and collagen reads and for elastin, reads with . Since ρ remains constant with representing the current mass density of constituent α within layer Γ (i.e., per unit current volume of the respective layer Γ at each time τ).

Passive and active stresses

The mechanical response of an artery is assumed to be isochoric for transient deformations at each fixed G&R time s. The layer-specific Cauchy stress tensor thus reads in the media, and in the adventitia, where is the deformation-dependent part of the Cauchy stress for constituent α within layer Γ, is the active stress tensor generated by the smooth muscle within the media, and pΓ are layer-specific pressure-type Lagrange multipliers associated with the incompressibility constraints and (transiently constant) that are to be determined from equilibrium and boundary conditions at fixed G&R times.

The passive Cauchy stresses for collagen and smooth muscle can be obtained from their associated second Piola–Kirchhoff stresses, deriving from the layer-specific strain energy functions of Eq. (13) with the second Piola–Kirchhoff stress tensor at the constituent level, deriving from , as The push-forward operation over yields the respective Cauchy stress tensor to be used further in Eqs. (17) or (18).

If we consider that elastin is neither produced nor degraded for s > 0 (because functional elastin is produced during the perinatal period, and elastin tends to degrade very slowly except in pathologies characterized by marked proteolytic activity) then and the strain energy function for elastin is The resulting second Piola–Kirchhoff stresses are where represents the layer-specific stress tensor at the constituent level. The passive Cauchy stresses in Eqs. (17) or (18) are obtained from Eq. (21), with α = e.

The active tensile stress generated by the smooth muscle tone within the media is considered to be exerted along the circumferential direction e, namely, where is the spatial mass fraction, is the maximum stress that the muscle can generate, C(s) > 0 is a ratio of vasoconstrictors (e.g., endothelin-1) to vasodilators (e.g., nitric oxide), λ and λ0 are the stretches at which the active force generating capability either is maximum or vanishes, respectively, and is the current active muscle fiber stretch. The ratio C(s) is written in terms of wall shear stress through where C is a basal ratio and C is a scaling factor, noting that vasodilators are produced by the endothelium when Δτ > 0 and vasoconstrictors are produced when Δτ < 0. Finally, the circumferential stretch for the active tone is defined as , with a(s) being the current luminal radius and a(s) being an active reference length whose evolution may be modeled through where 1/k is a characteristic time of remodeling and a(0) = a(0). An integral-type (convolution) solution of Eq. (26) for a(s) reads

Rate-independent solution for given pressure, flow rate, and axial stretch

We now recognize that the hereditary integral formulation accounts primarily for the finite half-lives of the cells and matrix that were incorporated within the evolving tissue at different past times. If, however, the tissue adapts fully to the prior perturbation in loading and all constituents produced during the adaptive period have been replaced with new constituents (continuously) produced in the new homeostatic state, then we can pre-integrate the prior, time-dependent G&R formulation following Ref. 13. In doing so, we obtain an associated time-independent solution for the bilayered G&R model outlined above. As we show below, this particularized formulation yields a system of algebraic equations that can be solved efficiently and whose solution represents either a steady-state, long-term solution reached at a G&R time s much greater than a characteristic time s over which the G&R processes take part, namely, at s ≫ s, or more importantly a quasi-equilibrium G&R evolution at any time s (in this case, s playing the role of a parameter), as we explain in Sec. IV C.

Substitution of Eq. (10) into Eq. (12) yields which can be integrated for constant (i.e., fully evolved homeostatic, denoted by subscript h) values at times s ≫ s, with and , to give since the integral of , as given in Eq. (7) with constant , is Hence, from Eq. (29), a mechanobiologically equilibrated G&R process associates with an equilibrium value of the mechano-stimulus function for collagen and smooth muscle in Eq. (10), namely, which, by virtue of Eq. (11), requires either Δσ = 0 = Δτ or, more generally As in Ref. 13, we take τ = 4μQ/(πa3), with Q being the volumetric flow rate and μ the blood viscosity, and in Eq. (9) as the first principal invariant of the mean wall Cauchy stress , namely, , where we assume a quasi-plane-stress state for which . The mean in-plane (biaxial) stresses σ and σ are given in terms of the distending pressure P and the global axial force on the vessel f, respectively, through The intramural and shear “over-stresses” expressed in terms of possibly fully evolved (h) or original (o) homeostatic values (i.e., we allow new homeostatic set-points to evolve) read where a (present in σ, σ and τ), h and h (in σ and σ), and f (in σ) are unknowns to be determined for each prescribed alteration in blood pressure, γ = P/P, blood flow, ε = Q/Q, and axial stretch λ = l/l (note that, in typical biaxial experiments, one usually prescribes axial stretch rather than axial load and it is not possible to infer f). For simplicity, let smooth muscle and collagen share the same ratio of gain parameters , despite generally different values of and . This means that the perturbation functions are proportional, hence the linearly dependent equations resulting from Eq. (32) reduce to

So far, we have a single equation with four unknowns, namely Eq. (35) in terms of a, h, h, and f. The layer-specific Jacobians J and J, expressed in terms of the “unchanging” homeostatic stretches in , introduce two additional unknowns (i.e., J and J), namely, where the radial and circumferential stretches are expressed in terms of a, h and h through Eqs. (3) and (4), and we assume a common prescribed axial stretch λ. Since the mass of elastin does not change in physiologic adaptations, its layer-specific spatial mass density becomes directly related to its original spatial mass density through the corresponding volume ratios which provide two additional equations but also two additional unknowns, namely and . Equation (16) particularized to both the media and adventitia yields two more equations, but introduces three more unknowns (namely , and ) which leaves four more unknowns than equations. However, as shown previously, additional relations in terms of the different (evolved homeostatic) spatial mass densities of collagen and smooth muscle can be written as and where , and . We also have, for different layer-specific cohorts i of collagen (e.g., circumferential, axial, and symmetric diagonal) where . Finally, the system of equations is closed by the global equilibrium equations σ = P and σ(2a + h) = f, where the internal circumferential force σ (per unit axial length) in terms of layer-specific stresses in Eqs. (17) and (18), yields and the internal axial force yields, similarly, The different expressions of the mechanobiologically equilibrated, layer-specific stresses to be used in Eqs. (17) and (18) are—see details in Ref. 13 and

Importantly, the nonlinear equations derived in this section do not depend on G&R time s, which means that they yield a mechano-adapted solution of the artery for a given sustained altered pressure P, flow rate Q, and axial stretch λ—hence constituting a truly time- and rate-independent G&R formulation for an idealized bilayered artery. In practice, this system of equations and unknowns may be reduced to five equations and unknowns as explained in Sec. IV E and illustrated in examples addressed above in Results (Sec. II). Of course, other resolution procedures are possible.

The quasi-equilibrium hypothesis

To arrive at the equilibrium condition of Eq. (31), from Eq. (28), we assumed that the artery preserves a static state for a sufficiently long time such that integration of Eq. (29) is exact. The same holds for the equilibrium stresses of cohorts of smooth muscle and collagen given in Eq. (44), which are obtained upon integration of the corresponding integral-type Eqs. (19) and (21). Considering Eq. (29), we see that and , as given in Eqs. (8) and (11), become constant if the wall stress metric and the shear stress τ remain constant over long times. Taking into account Eqs. (33) and (34), this requires that the applied external loads P, Q, and λ remain constant over long times, which are, indeed, the ultimate variables that stimulate the G&R response in the present mechanoadaptive case.

In what follows, we relax this equilibrium hypothesis to introduce the definition of a quasi-equilibrium state, which we illustrate by comparing three cases: I, II, and III. Assume a survival function , as given in Eq. (7), as a function of the deposition time τ ≤ s0 for a fixed G&R time s0. Let be a proper integration domain such that , with 0+ a sufficiently small value such that longer integrations of beyond s0–Δs yield (additional) negligible contributions. Consider, first, the evolution over time of a generic nondimensional external insult that, qualitatively, includes any mechanical alteration such as , or . Obviously, if the rate of change of and are comparable within , namely, for , then the associated variations of and over time [assuming gain parameters and of order unity, see Eq. (11)] cannot be neglected in the integral of Eq. (28), and the general formulation cannot be simplified. In this case (I), one needs to keep track the history of all the involved variables, at least within , to compute the G&R solution at any time s0—that is, one must employ the full (hereditary integral) constrained mixture formulation.

Consider here, however, a case (II) for which within . Then, if the G&R response is mechanobiologically stable and has relaxed fully, then the variables and reach, after the corresponding characteristic G&R period, constant values and within the relevant integration domain as well. The artery thus reaches a (new) tissue maintenance state, for which is also constant, leading to the result of Eq. (29) and, eventually, to the mechanobiologically equilibrated formulation detailed above. As explained in, a characteristic time scale for both mass removal and production processes is given by , hence mechanobiological equilibrium is attained at s ≫ s. Recall that gain parameters and can modify the value of s quantitatively.

Another important case (III) is one for which the rate of change of the external stimuli is much less than the rate of change of the removal function within . From a mathematical standpoint, the variables of any integrands can then be regarded constant within the integration domain , hence we recover the previous equilibrium formulation approximately at each G&R time s, namely a quasi-equilibrium formulation valid during the evolution. For example, the integral of Eq. (28) yields, for slow , at any G&R time s [cf. Eq. (29)] whereby the stimulus function resolves, that is, [cf. Eq. (31)] From a physical standpoint, we can equivalently focus on the evolution of an arbitrary mass deposited at a fixed time τ0 that is removed gradually during . Hence, if quasi-equilibrium conditions are satisfied, the differential mass senses an almost constant mechanical environment during its timespan Δs, undergoing then a so-called quasi-steady evolution. The evolution is said to be mechanobiologically quasi-equilibrated (rather than mechanobiologically equilibrated) because the external stimuli may be different at different G&R times over a time scale much longer than Δs; that is, for , but for s1 – s0 ≫ Δs, in general.

In summary, if within , the time-independent formulation derived in Sec. IV B remains valid throughout the mechanoadaptation, with G&R time s playing the role of a parameter while the external stimuli is able to change over a much longer period. Using the decay function of Eq. (7), we obtain (by the Leibniz integral rule) so a characteristic mass-specific rate of removal at a given G&R time s is . If we include the evolution of the active reference length for the active stress contribution of smooth muscle, the quasi-equilibrium formulation is valid then if . Note that this is equivalent to saying that a characteristic time for G&R, namely s = 1/k, is much shorter than a characteristic time of change of the (normalized) external loads, namely s = 1/k.

Temporal nondimensionalization of the full model

The concept of a characteristic time scale for a G&R response of a living soft tissue encourages us to nondimensionalize in the time domain in the full constrained mixture model outlined in Sec. IV A. Given the characteristic time , we first define dimensionless rate parameters and as well as dimensionless time variables such that we have equivalent products between original and dimensionless rate parameters and time variables of the type and so forth. We can then obtain equivalent integrals for Eq. (7) where [cf. Eq. (8)] and for Eq. (12) where [cf. Eq. (10)] and for Eq. (13) Finally, the time-dimensionless counterpart of the evolution equation in the rate form of Eq. (26) reads Thus, every time-dependent evolution equation can be equivalently expressed in terms of dimensionless rate parameters and time variables. This means that, once a time-dependent G&R response is computed for a given tissue, say tissue B, with characteristic G&R time , rate-type parameters k and , and prescribed external insults , the solution represented in terms of the dimensionless time is common for any other tissue, say D, with proportional rate-type parameters and , with , and prescribed external insults , with common remaining material properties, because they share the same time-dimensionless formulation.

Resolution procedure

For the pre-integrated model of Sec. IV B, we will solve the system of nonlinear equations formed by Eqs. (35), (38)1, (38)2, (42), and (43), where the unknowns are the (potentially new) homeostatic inner radius a, layer thicknesses h and h, spatial mass density of collagen within the media , and global axial force f, with other variables expressed easily in terms of the selected unknowns. The resulting system of equations is time-independent and so too its outcome.

The hyperelastic mechanical response of elastin is modelled using a neoHookean relation with c being the shear modulus. Hyperelastic responses of both smooth muscle and collagen are modelled using Fung-type relations where (dimensions of stress) and (dimensionless) are material parameters, and is the corresponding fiber stretch. We consider four collagen fiber families in both the media and adventitia: one oriented circumferentially (labelled with θ), one oriented axially (z), and two oriented in symmetric diagonal (d) directions  ± α0 with respect to the axial direction. The contributions of circumferential collagen and smooth muscle are combined in the media (referred to as medial circumferential smooth muscle m). Medial and adventitial collagen are assumed to share the same hyperelastic, rate, and gain constants, the difference in contributions coming from different mass fractions. All the material parameters needed to obtain both time-dependent and time-independent solutions are listed in Table I. The specific values of the parameters are best-fit values determined from in vitro biaxial data from passive elastic arteries of mice. In order to show the full consistency between both formulations when they include all possible contributions to stress, however, we take additional values for the active response of smooth muscle (Table I). Finally, no ethics approval was required since all work was numerical.