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.


I. INTRODUCTION
As aptly stated in 1995 by Fung, "Every specialty in biomechanics begins with the study of constitutive relations," 1 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. 2 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 a) Electronic mail: m.latorre.ferrus@upm.es b) Author to whom correspondence should be addressed: jay.humphrey@yale.edu 2473-2877/2018/2(2)/026108/20 V C Author(s) 2018. 2, 026108-1 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, mechanoadaptation 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, [3][4][5][6][7][8] as well as other cellular and tissue-level processes. [9][10][11] 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. 2 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 12 while another has been to derive an associated steady-state form that reveals the final (evolved) state, 13 not unlike equivalence derivations in viscoelasticity that relate integral and rate forms 14 or those that focus on long-term responses in relaxation and creep. 15 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 rateindependent ("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.

II. RESULTS
A. 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 thinwalled 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).
Panel (b) shows the resulting/evolving stimulus function ! c 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 q c CR , a, h M , and h A , 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 ! a C ) 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 G&R ¼ 1/7 % 0.143 days À1 ) is greater than the characteristic rate of change of the external loads (in this example, k ext $ 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 ! c [ Fig. 1(b)] and smooth muscle ! m (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 G&R ) s G&R ¼ 7 days, when mechanobiological equilibrium is restored [mathematically described by ! a h ¼ 1, 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 (s wh 6 ¼ s wo , r hhh 6 ¼ r hho , and r zzh 6 ¼ r zzo ) at the new evolved state, which, yet satisfy the more general equilibrium condition ! a h ¼ ! a o ¼ 1. 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 c ¼ P h /P o ¼ e ¼ Q h /Q o ¼ 1.15, here computed via a Newton-Raphson method with initial guess given by an ideal adaptation of the type This solution is shown by solid squares in Fig. 1, which reveals precise correspondence with the full, timedependent constrained mixture model at long (fully adapted) times.
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 timedelayed 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.

B. 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 G&R is much greater than the characteristic rate of change of the external stimuli k ext , satisfying then the quasi-steady-state relation k ext / k G&R ( 1 at any G&R time s. For these particular situations, the arterial adaptation (occurring over a time scale s G&R ¼ 1/k G&R ) to different alterations in mechanical loading (occurring over a time scale s ext ¼ 1/k ext ) 1/k G&R ¼ s G&R ) 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 G&R should be with respect to the external loading characteristic time s ext 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 G&R ¼ 1/k G&R ¼ 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 c(s) P(s)/P o shown in Fig. 2(a). In order to assess the quasisteady-state assumption for each case, these values of s G&R 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 o ! 1.15, namely s ext $ 10 days in this case. We also show in the same figure the single rateindependent solution of Sec. IV B (pseudoelastic, solid line) computed as a function of the time-dependent inner pressure ratio P(s)/P o where time s simply plays the role of a simulationdriver parameter. In this case, we start the iterative solution procedure at each new time step using the converged solution at the previous time step. Figure 2 shows that the prescribed increase in pressure provokes simultaneous changes in referential mass densities q c MR and q c AR  Fig. 2(b) for collagen a ¼ c, with similar tendency for smooth muscle], in four simulations (three rate-dependent based on different characteristic times s G&R and one rate-independent). In particular, the evolution of inner radius, thicknesses, and masses 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 G&R ¼ {0.1, 1, 10} days, for an isolated increase in pressure with s ext $ 10 days. Shown are (a) prescribed load P/P o from 1 to 1.15, (b) mechano-stimulus function ! c , (c) referential mass densities of collagen q c MR =q c Mo ¼ q c AR =q c Ao , (d) relative inner radius a/ a o , (e) relative medial thickness h M /h Mo , and (f) relative adventitial thickness h A /h Ao . 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 longterm 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. for s G&R /s ext $ 0.1/10 ¼ 0.01, for which ! a ' 1 [cf. Eq. (48)], is (in practice) indistinguishable from the rate-independent evolution, for which ! a h ¼ 1 [cf. Eq. (31)]. For the other two simulations, s G&R /s ext $ 1/10 ¼ 0.1 and s G&R /s ext $ 10/10 ¼ 1, the mechano-stimulus function no longer satisfies the quasi-equilibrium condition ! a ' 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, 22 i.e., a/a o ' 1 ¼ e 1=3 where e ¼ Q/Q o . Note, however, that the overall solution (including medial and adventitial thicknesses and constituent masses) given by the rateindependent 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 o ¼ 1.23 which is 7% greater than that for an ideal mechanoadaptation 1.15 ¼ ce 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. 23 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 G&R /s ext $ 1). Indeed, the increase in pressure mainly provokes an initial increase in circumferential stress, hence the stimulus functions ! a drive a growth process (with parallel increments of thickness and mass) until, eventually, stresses return close to normal values and ! a ! ! a h ¼ 1. Figure 3 shows the same type of analysis but for a particular increase in flow rate e(s) . 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 s G&R =s ext $ 0:01; 8s) and exactly the same long-term, tissue-maintenance solution as the full model in all cases (s G&R =s ext $ 0:01; 0:1; 1 f g , s ¼ 100 days). If we compare the full model predictions for s G&R /s ext $ 0.1 and s G&R /s ext $ 1 to the single one given by the rateindependent 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, a=a o ¼ 1:04 ' 1:05 ¼ e 1=3 and h=h o ' 1:01 < 1:05 ¼ ce 1=3 (not shown), indicative of a near but not fully mechanoadaptive solution again. In the case of comparable timescales (s G&R /s ext $ 1), the increase in flow rate mainly provokes an increase in shear stress, and thus a decrease in the stimulus functions ! a [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 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). consistent with vasodilation with nitric oxide), provokes an increase in circumferential stress that, subsequently, drives a growth process via ! a > 1. Finally, stresses closely return to normal values, such that ! a ! ! a h ¼ 1. Similar conclusions regarding quasi-equilibrium (s G&R =s ext $ 0:01; 8s) and long-term (s G&R =s ext $ 0:01; 0:1; 1 f g , s ¼ 100 days) predictions are obtained for a single increase in the axial stretch ratio k z (s)/k zo (in this case up to k zh /k zo ¼ 1.10, Fig. 4). Regarding short-and mid-term evolutions, especially for s G&R /s ext $ 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. 24,25 Finally, the common mechanobiologically equilibrated state (not fully reached at s ¼ 100 days for s G&R /s ext $ 1) is such that a/a o ' 1 ¼ e 1=3 and h/h o ¼ 1.05 > 1 ¼ ce 1=3 . In this case, for comparable timescales (s G&R /s ext $ 1), the increase in axial stretch mainly provokes an increase in axial stress and thus stimulus functions ! a , 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 ! a < 1. Finally, because of the decreased thickness, circumferential stress increases slightly, driving G&R until the artery reaches a mechanobiologically equilibrated state associated with ! a ! ! a h ¼ 1. Interestingly, the mechano-stimulus functions ! a [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 c(s) ¼ P(s)/P o for the slowest response examined in Fig. 2 (s G&R / s ext $ 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.

III. 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 FIG. 4. Similar to Fig. 2 except for an isolated increase in axial stretch k z 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). 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. [3][4][5][6][7][8][9] 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 13 and quasiequilibrium 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 G&R , 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. 15,[26][27][28] In parallel, time-independent G&R models as presented herein will never replace time-dependent G&R models, either fully integrated 2 or kinematically motivated. 29 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 2 and mechanobiologically equilibrated 13 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 G&R $ 0.1 days), which we used as a more stringent test of the concept. Nevertheless, values of s G&R should be assessed relative to a characteristic time of changes in the loads that stimulate G&R, namely s ext . Thus, a dimensionless number of the type s G&R /s ext , comparing characteristic times of G&R and loading, is what actually determines the goodness of a quasi-static assumption. Regarding actual arterial behaviors, s G&R may range from $10 to 100 days 17,33 whereas changes in loading may range from seconds to many days, 21,34 hence dimensionless numbers s G&R /s ext of order unity, or lower, exist for s ext $10 to 100 days or longer. The lower the value of s G&R /s ext , the better the pseudoelastic G&R approximation, as one would expect from viscoelasticity theory. In contrast, the greater the ratio s G&R /s ext , the better the transient "elastic" (without G&R) approximation; see Example 2 in Ref. 13. Because s G&R is material-dependent, it can change during certain conditions or diseases, 33 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. 35 Indeed, kinematic models of growth 29 ultimately rely on the recognition of these two time scales, with evolution of the growth tensor F g depending on the "growth (or) remodeling timescale." 35 Hence, growth laws may be expressed in kinematic models in terms of rate parameters 36 that, alternatively, may depend on problemspecific growth metrics to prevent unlimited growth, 37 with baseline values identified as rates of initial growth. 38 Similarly, the so-called global growth approach 39 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.

A. A bilayered constrained mixture model for arteries
First consider the kinematics and evolution equations of a constrained mixture model 2 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. 1. 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. 16 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 j o ¼ j(0), serve as a reference. The configuration at the current G&R time s > 0 is denoted j(s), while configurations at all intermediate times 0 s s are denoted by j(s). Hereafter, for the sake of notational simplicity, we will omit the time dependence when it is not needed explicitly.  17 but we consider here a faster rate of turnover 33 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, h, 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 r and s denote intramural and wall shear stresses, respectively. See Table II   , showing the different deformation gradients F M and F A experienced by the medial and adventitial layers over time, in particular at times 0 s s. Shown, also, are the different natural configurations j a nðsÞ for the different structurally significant constituents a (elastin "e," smooth muscle "m," and collagen "c"), which are deposited with separate but constant deposition stretches G a within both layers at the indicated deposition times s (except for the smooth muscle, present in media only), as well as geometric parameters a (luminal radius), h M (medial thickness), and h A (adventitial thickness) and luminal pressure P, which can also change over time.
We denote the inner (luminal) radius of the artery in a generic configuration as a, the medial thickness as h M , and the adventitial thickness as h A . The total wall thickness is thus h ¼ h M þ h A . 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 MA ¼ a þ h M , hence l M ¼ l A l.
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 j o and a generic configuration j. Using cylindrical coordinates X rhz ¼ {r, h, 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 MA ¼ a þ h M 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 / M ¼ ½/ e M ; / m M ; / c M satisfying P e;m;c a / a M ¼ 1 in the media, and / A ¼ ½/ e A ; / m A ; / c A satisfying P e;m;c a / a A ¼ 1 in the adventitia, with / m A ¼ 0. 17 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 M and F A in Eqs. (1) and (2). Each constituent and cohort, however, has its own evolving natural configuration j a n ðsÞ, where s 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 a , which means that the natural configuration of a constituent a at time s, namely, j a n ðsÞ, is common for both layers (see Fig. 6). In addition, we assume (in physiologic adaptations) constant and volume-preserving deposition stretch tensors G a , with detG a ¼ 1. Following the deformation path shown in Fig. 6, the deformation gradient F a CnðsÞ ðsÞ experienced by constituent a deposited at time s within layer C ¼ 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 4 and for elastin, which is produced perinatally and thus at s 0 in our case of G&R in maturity, 6 as since F C (s ¼ 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 q a C ðs; sÞ 2 ½0; 1 of the type (we use the form introduced in Ref. 13, but others are possible 18 ) where the rate "parameter" k a C ðtÞ, with s t s, is assumed (constitutively) to increase with respect to its original homeostatic constant value k a Co through with Dr(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, namelyrðtÞ, and its corresponding homeostatic valuer o , 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 a per unit reference volume of the mixture (in this case, within layer C), we have for both layers (see Tables I and II for where m a CN ðsÞ ¼ k a C ðsÞq a CR ðsÞ > 0 is an evolving nominal mass production rate including, importantly, the same function k a C ðsÞ employed for mass removal and the referential mass density q a CR ðsÞ of constituent a within layer C (i.e., per unit reference volume of the respective layer C). This relation ensures balanced production and removal in homeostatic states, 13 for which m a CR ! m a CN . Moreover, ! a C ðsÞ 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 ! a C ðsÞ for arteries subjected to perturbed blood pressures and flows, linearized about the initial homeostatic state, reads 4 where K a Cr and K a Cs are constituent-and layer-specific (constant) gain parameters, s w is the flow-induced shear stress acting over the endothelium, and Ds w ¼ (s w -s wo )/s wo . In particular, ! a Co ¼ 1 at the original homeostatic state, where Dr ¼ 0 ¼ Ds w . Given constitutive equations for production and removal, the evolution of the layer-specific referential mass density of the cohort a is 4,13 Assuming that the spatial mass density q of the arterial wall (i.e., its current mass per unit current volume of mixture) remains constant, the strain energy function of constituent a within layer C, also defined per unit reference volume of layer C, reads 13 whereŴ a ðC a CnðsÞ ðsÞÞ is the volume-specific strain energy function of constituent a and C a CnðsÞ ðsÞ is the right Cauchy-Green deformation tensor obtained from F a CnðsÞ ðsÞ, which for smooth muscle and collagen reads and for elastin, reads with with q a C representing the current mass density of constituent a within layer C (i.e., per unit current volume of the respective layer C at each time s).

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 r a C is the deformation-dependent part 5 of the Cauchy stress for constituent a within layer C, r act is the active stress tensor generated by the smooth muscle within the media, and p C are layer-specific pressure-type Lagrange multipliers associated with the incompressibility constraints J M ¼ detðF M Þ and J A ¼ detðF A Þ (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 W a CR of Eq. (13) 13 The push-forward operation over S a C s ð Þ yields the respective Cauchy stress tensor r a C ðsÞ 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 19 ) then q e CR ðsÞ ¼ q e CR ð0Þ ¼: q e CR and the strain energy function for elastin is The resulting second Piola-Kirchhoff stresses are whereŜ e C ¼ 2@Ŵ e ðC e C Þ=@C e C 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 a ¼ e.
The active tensile stress generated by the smooth muscle tone within the media is considered to be exerted along the circumferential direction e h , namely, 3 where / m M ðsÞ ¼ q m M ðsÞ=q is the spatial mass fraction, T max 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), k M and k 0 are the stretches at which the active force generating capability either is maximum or vanishes, respectively, and k mðactÞ h ðsÞ is the current active muscle fiber stretch. The ratio C(s) is written in terms of wall shear stress through 3 C s ð Þ ¼ C B À C S Ds w s ð Þ; where C B is a basal ratio and C S is a scaling factor, noting that vasodilators are produced by the endothelium when Ds w > 0 and vasoconstrictors are produced when Ds w < 0. Finally, the circumferential stretch for the active tone is defined as k mðactÞ h ðsÞ ¼ aðsÞ=a act ðsÞ, with a(s) being the current luminal radius and a act (s) being an active reference length whose evolution may be modeled through 3 da act s ð Þ ds ¼ k act a s ð Þ À a act s ð Þ À Á ; where 1/k act is a characteristic time of remodeling and a act (0) ¼ a(0). An integral-type (convolution) solution of Eq. (26) for a act (s) reads k act a s ð Þe Àk act ðsÀ sÞ d s : B. 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, timedependent G&R formulation following Ref. 13. In doing so, we obtain an associated timeindependent 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 G&R over which the G&R processes take part, 13 namely, at s ) s G&R , 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 G&R , with k a C ! k a Ch ; q a CR ! q a CRh and ! a C ! ! a Ch , to give since the integral of q a C ðs; sÞ ! q a Ch ðs; sÞ, as given in Eq.
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 Dr h ¼ 0 ¼ Ds wh or, more generally K a Cr Dr h À K a Cs Ds wh ¼ 0 : As in Ref. 13, we take s w ¼ 4lQ/(pa 3 ), with Q being the volumetric flow rate and l the blood viscosity, andr in Eq. (9) as the first principal invariant of the mean wall Cauchy stress r, namely,r ¼ tr r ' r hh þ r zz , where we assume a quasi-plane-stress state for which r rr =r ' 0. The mean in-plane (biaxial) stresses r hh and r zz are given in terms of the distending pressure P and the global axial force on the vessel f z , 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 Dr h ¼ r hhh þ r zzh r hho þ r zzo À 1; and Ds wh ¼ s wh s wo À 1; where a h (present in r hhh , r zzh and s wh ), h Mh and h Ah (in r hhh and r zzh ), and f zh (in r zzh ) are unknowns to be determined for each prescribed alteration in blood pressure, c h ¼ P h /P o , blood flow, e h ¼ Q h /Q o , and axial stretch k zh ¼ l h /l o (note that, in typical biaxial experiments, one usually prescribes axial stretch rather than axial load and it is not possible to infer f z in vivo 17 ). For simplicity, let smooth muscle and collagen share the same ratio of gain parameters Cs , despite generally different values of K a Cr and K a Cs . This means that the perturbation functions ! a C ðsÞ À 1 are proportional, 13 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 Mh , h Ah , and f zh . The layer-specific Jacobians J Mh and J Ah , expressed in terms of the "unchanging" homeostatic stretches in ½F Ch rhz ¼ diag ½k Crh ; k Chh ; k Czh , introduce two additional unknowns (i.e., J Mh and J Ah ), namely, where the radial and circumferential stretches are expressed in terms of a h , h Mh and h Ah through Eqs. (3) and (4), and we assume a common prescribed axial stretch k zh . Since the mass of elastin does not change in physiologic adaptations, its layer-specific spatial mass density q e Ch becomes directly related to its original spatial mass density q e CR q e Co through the corresponding volume ratios which provide two additional equations but also two additional unknowns, namely q e Mh and q e Ah . Equation (16) particularized to both the media and adventitia yields two more equations, but introduces three more unknowns (namely q m Mh ; q c Mh , and q c Ah ) q e Mh þ q m Mh þ q c Mh ¼ q; and q e Ah þ q c Ah ¼ q; which leaves four more unknowns than equations. However, as shown previously, 13 additional relations in terms of the different (evolved homeostatic) spatial mass densities of collagen and smooth muscle can be written as (40) Ms . We also have, for different layer-specific cohorts i of collagen (e.g., circumferential, axial, and symmetric diagonal) 13 where q c C ¼ P q c i C . Finally, the system of equations is closed by the global equilibrium equations r hhh h h ¼ P h a h and r zzh ph h (2a h þ h h ) ¼ f zh , where the internal circumferential force r hhh h h (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 h , flow rate Q h , and axial stretch k zh -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.

C. 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 k a C and ! a C , as given in Eqs. (8) and (11), become constant if the wall stress metricr r hh þ r zz and the shear stress s w remain constant over long times. Taking into account Eqs. (33) and (34), this requires that the applied external loads P, Q, and k z 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 quasiequilibrium state, which we illustrate by comparing three cases: I, II, and III. Assume a survival function q a 0 ðsÞ ¼ q a C ðs 0 ; sÞ, as given in Eq. (7), as a function of the deposition time s s 0 for a fixed G&R time s 0 . Let T a :¼ ½s 0 À Ds a ; s 0 be a proper integration domain such that q a 0 2 ½0 þ ; 1, with 0 þ a sufficiently small value such that longer integrations of q a 0 ðsÞ beyond s 0 -Ds a yield (additional) negligible contributions. Consider, first, the evolution over time of a generic nondimensional external insult n I ext ðsÞ that, qualitatively, includes any mechanical alteration such as PðsÞ=Pðs 0 Þ 6 ¼ 1; QðsÞ=Qðs 0 Þ 6 ¼ 1, or k z ðsÞ=k z ðs 0 Þ 6 ¼ 1. Obviously, if the rate of change of q a 0 ðsÞ and n I ext ðsÞ are comparable within T a , namely, _ q a 0 ¼ dq a 0 =ds $ dn I ext =ds ¼ _ n I ext for s 2 T a , then the associated variations of k a C ðsÞ and ! a C ðsÞ over time [assuming gain parameters K a Cr and K a Cs 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 T a , to compute the G&R solution at any time s 0 -that is, one must employ the full (hereditary integral) constrained mixture formulation. 2,18 Consider here, however, a case (II) for which _ n II ext ¼ 0 within T a . Then, if the G&R response is mechanobiologically stable 20 and has relaxed fully, then the variables k a C and ! a C reach, after the corresponding characteristic G&R period, constant values k a Ch and ! a Ch ¼ 1 within the relevant integration domain T a as well. The artery thus reaches a (new) tissue maintenance state, for which q a CR ! q a CRh is also constant, leading to the result of Eq. (29) and, eventually, to the mechanobiologically equilibrated formulation detailed above. As explained in, 13 a characteristic time scale for both mass removal and production processes is given by s G&R ¼ 1=minfk a Co ; k act g, hence mechanobiological equilibrium is attained at s ) s G&R . Recall that gain parameters K a Cr ¿1 and K a Cs ¿1 can modify the value of s G&R quantitatively. Another important case (III) is one for which the rate of change of the external stimuli _ n III ext is much less than the rate of change of the removal function _ q a 0 within T a . From a mathematical standpoint, the variables of any integrands can then be regarded constant within the integration domain T a , 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) From a physical standpoint, we can equivalently focus on the evolution of an arbitrary mass deposited at a fixed time s 0 that is removed gradually during T a :¼ ½s 0 ; s 0 þ Ds a . Hence, if quasi-equilibrium conditions are satisfied, the differential mass senses an almost constant mechanical environment during its timespan Ds a , 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 n III ext may be different at different G&R times over a time scale much longer than Ds a ; that is, n III ext ðs 0 Þ ' n III ext ðs 1 Þ for s 1 À s 0 $ Ds a , but n III ext ðs 0 Þ6 'n III ext ðs 1 Þ for s 1s 0 ) Ds a , in general. In summary, if _ n ext ( _ q a 0 within T a , 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 _ q a 0 ðsÞ ¼ k a C ðsÞ $ k a Co . 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 _ n ext ðsÞ k ext ðsÞ ( k G&R ¼ minfk a Co ; k act g; 8s. Note that this is equivalent to saying that a characteristic time for G&R, namely s G&R ¼ 1/k G&R , is much shorter than a characteristic time of change of the (normalized) external loads, namely s ext ¼ 1/k ext .

D. 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 s G&R ¼ 1=minfk a Co ; k act g, we first define dimensionless rate parametersk a Co ¼ k a Co Á s G&R ; a ¼ m; c; 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 s B G&R , ratetype parameters k actB and k aB Co ; a ¼ m; c, and prescribed external insults n B ext ðs B Þ, the solution represented in terms of the dimensionless times ¼ s B =s B G&R is common for any other tissue, say D, with proportional rate-type parameters k actD ¼ bk actB and k aD Co ¼ bk aB Co , with b ¼ s B G&R =s D G&R , and prescribed external insults n D ext ðs D Þ ¼ n B ext ðs B =bÞ, with common remaining material properties, because they share the same time-dimensionless formulation.

E. Resolution procedure
For the pre-integrated model of Sec. IV B, we will solve the system of nonlinear equations formed by Eqs.
where c a 1 (dimensions of stress) and c a 2 (dimensionless) are material parameters, and k a nðsÞ ðsÞ is the corresponding fiber stretch. We consider four collagen fiber families in both the media and adventitia: one oriented circumferentially (labelled with h), one oriented axially (z), and two oriented in symmetric diagonal (d) directions 6 a 0 with respect to the axial direction. The contributions of circumferential collagen and smooth muscle are combined in the media 16 (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. 21 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.