Aortic Remodeling Kinetics in Response to Coarctation-Induced Mechanical Perturbations

Abstract Background: Coarctation of the aorta (CoA; constriction of the proximal descending thoracic aorta) is among the most common congenital cardiovascular defects. Coarctation-induced mechanical perturbations trigger a cycle of mechano-transduction events leading to irreversible precursors of hypertension including arterial thickening, stiffening, and vasoactive dysfunction in proximal conduit arteries. This study sought to identify kinetics of the stress-mediated compensatory response leading to these alterations using a preclinical rabbit model of CoA. Methods: A prior growth and remodeling (G&R) framework was reformulated and fit to empirical measurements from CoA rabbits classified into one control and nine CoA groups of various severities and durations (n = 63, 5–11/group). Empirical measurements included Doppler ultrasound imaging, uniaxial extension testing, catheter-based blood pressure, and wire myography, yielding the time evolution of arterial thickening, stiffening, and vasoactive dysfunction required to fit G&R constitutive parameters. Results: Excellent agreement was observed between model predictions and observed patterns of arterial thickening, stiffening, and dysfunction among all CoA groups. For example, predicted vascular impairment was not significantly different from empirical observations via wire myography (p-value > 0.13). Specifically, 48% and 45% impairment was observed in smooth muscle contraction and endothelial-dependent relaxation, respectively, which were accurately predicted using the G&R model. Conclusions: The resulting G&R model, for the first time, allows for prediction of hypertension precursors at neonatal ages that is currently challenging to examine in preclinical models. These findings provide a validated computational tool for prediction of persistent arterial dysfunction and identification of revised severity–duration thresholds that may ultimately avoid hypertension from CoA.


Introduction
Coarctation of the aorta (CoA) is a constriction of the proximal descending thoracic aorta and is one of the most common congenital cardiovascular defects [1]. Surgical repair remains the standard approach for treatment in most young infants, while catheter intervention with balloon angioplasty and/or stenting is the standard of care in older children and adults [2,3]. Guidelines for intervention include a transcatheter peak-to-peak blood pressure (BP) gradient ≥20 mmHg, and many published reports regard this criterion as a hemodynamically significant CoA [3]. Unfortunately, hypertension and left-ventricular hypertrophy are common despite successful CoA repair based upon this threshold [4,5].
Mechanistic understanding of hypertension in CoA is hindered by the paucity of data to investigate changes in the kinetics of growth and remodeling (G&R) associated with coarctation-induced mechanical stimuli [6,7]. Hence, animal models are used to investigate underling mechanisms and precursors of hypertension in response to CoA [8].
In contrast to the current putative threshold, recent studies using a rabbit model of CoA suggest that precursors of hypertension are present among treated rabbits with CoA gradients <20 mmHg [8][9][10]. This finding raises doubts about the efficacy of the current guidelines in preventing refractory hypertension after treatment. Interestingly, the duration of coarctation-induced mechanical stimuli also correlates with hypertension precursors [7,9], further supporting the importance of earlier treatment that has been shown to reduce hypertension prevalence among treated CoA patients [11].
In general, the arterial wall is assumed to adapt to mechanical stimuli (perturbation from normal state) and recover a homeostatic state via smooth muscle (SM) cell synthesis and matrix turnover [12][13][14]. For example, in CoA, the presence of the narrowing can cause elevated BP proximally and a turbulent velocity jet distally. This perturbation in blood pressure-induced intramural wall stress (IWS: σ) proximally and blood flow-induced wall shear stress (WSS: φ) distally changes the vascular tone toward recovering IWS and WSS back to preferred ranges that are assumed to be 'known' to the local vasculature. Consequently, fibrillar proteins such as collagen and SM are deposited to remodel the extracellular matrix and modify thickness, arterial caliber, and material properties. The kinematics of this compensatory thickening/stiffening response is stress-driven and has been previously quantified as a linear function of deviations from homeostatic values [14][15][16].
The aim of this work is to quantify G&R parameters applicable to CoA using a computational approach tuned to replicate arterial thickening, stiffening, and dysfunction measured in rabbit aortas exposed to mechanical stimuli within ranges seen clinically. Completion of the study objective will provide a validated computational tool for prediction of persistent arterial dysfunction and allow for the much-needed identification of severity-duration thresholds that prevent arterial dysfunctions contributing to refractory hypertension among treated CoA patients.

Stress-Mediated G&R Formulation
For a section of the aorta, at any time t, the current configuration can be represented simply as a pressurized cylindrical thin-wall tube. Deformation can, therefore, be described by mapping the reference point (R, Θ, Z) onto a new position (r, θ, z) in the same coordinate system when deformed: for some function r(R) and parameter κ. That is, deformations are restricted so that the principal axes of deformation coincide with the coordinate directions of the principal stretches.
Equilibrium in the radial direction gives the following [17]: with σ representing the IWS in the corresponding directions. Thickness integration with the σ r |r i = p as a boundary condition on the luminal surface (r = r i ), gives This equation can be rewritten in reference configuration by observing that r = λR and dr = (λκ) −1 dR from Equations (2) and (3), with the constitutive relation [17,18] σ θ − σ r = λ ∂W/∂λ: Neglecting through thickness variation of the integrand with the thin-wall assumption, i.e., h(t) << r i (t), yields where W(t) is the total strain energy per unit area at time t representing passive forces in the vasculature. On the other hand, W a (t) is the SM work density representing active forces acting on the arterial wall. Strain energy per unit area can be described by modeling the artery's ability to adapt to changes of mechanical environment to maintain a homeostatic state via vascular G&R [16]. This process occurs through removal of old vascular constituents and incorporation of new constituents, which can be described in the following time evolution form: where M k (t) describes the current mass per unit area for the vascular constituent k. The first term on the right side of the equation is the natural turnover of the initial mass. Q k (t) is the current remaining fraction for constituent k produced at time 0. The second term describes the natural turnover of the newly produced constituents, with q k (t − τ) being the remaining fraction of the constituent k produced at time τ and mass production rate m k (τ). Here, it is assumed that Q k (t) and q k (t) follow an exponential decay, i.e., exp −ν k t with rate constant ν k capturing turnover of constituent k, i.e., 1/80[days] for collagen and SM [19]. Elastin is assumed to neither degrade nor grow [20], i.e., Q e (t) = 1 and q e (t) = 0. For a constrained mixture, based on the mass-average principle, the total strain energy per unit area at time t is where W k (t, 0) and W k (t, τ) denote the strain energy density that constituent k ∈ (c, m, e) produces at times 0 and τ, respectively. The vessel is assumed a constrained mixture of collagen, SM, and elastin, i.e., c, m, e. Active stress in the circumferential direction σ a (t) can be defined from a constitutive relation assuming that vascular SM is oriented nearly circumferentially [21] with zero active stress in the radial and longitudinal directions.
The active contribution of circumferential stress is modulated as a linear function of generated stress per unit relaxed SM area, i.e., σ a (t) = λ m S. Here, G&R takes place on a relatively long timescale where the SM response is more importantly characterized by the phenotypic dedifferentiation between a contractile and a synthetic state [9,22,23] as compared to short-term calcium-driven cyclic interactions. This phenotypic shift has been reported as increased non-muscle myosin expression with decreased SM myosin [9] in response to coarctation-induced mechanical stimuli, including those derived from increased systolic, mean, and pulse BP in regions proximal to the coarctation. Although these arterial alterations are in the direction of increased medial thickness through sustained remodeling, the active contractile force is somewhat impaired possibly due to a phenotypic shift from the contractile to synthetic state [9]. Moreover, arterial thickening is at the expense of increased structural stiffness leading to augmented forward and backward pressure waves that further increase pulse pressure and BP [24,25]. This creates a positive feedback loop that may become irreversible leading to refractory hypertension even after removal of the CoA.
We assume herein that phenotypic dedifferentiation is a continuous process in which the SM cell gradually gains or loses its contractile response according to the function S = S (C), where C ≥ 0 is the phenotypic modulation stimulus. It represents the accumulation of stimuli descending from alterations in the stretch field. The generated active stress S is, therefore, associated with the slowly varying tone under normal physiological conditions [26]. Since not all infinitesimal areas of SM experience the same stretch, the expression for active stress is adapted to a situation with variable λ m , so that the total active stress difference σ a (t) − σ a r (t) is computed as a sum of active stress difference contributions from remaining infinitesimal area fractions produced at different times. To derive an equation for active stress differences, the derivative of Equation (9) is taken with respect to λ for SM constituent only. Multiplying both sides by λ then gives where the notation d f (x) = d f /dx. Substituting Equation (10) gives where λ m (t, 0) and λ m (t, τ) are the stretch ratio for SM constituents generated at times 0 and τ, respectively, and can be described in terms of stretch developed during growth, i.e., λ(t)/λ(τ), and pre-stretch of SMĜ m that develops at the time of production, which is assumed to be a constant equal to 1.3 according to Beak et al. [27,28]: The phenotypic modulation of contractile activity is taken to be Biomedicines 2023, 11, 1817

of 22
where S M is the maximum stress capacity of SM in its fully contractile state and assumed to be a constant equal to 100 kPa [29]. For a target lumen radiusr i , the stretch isλ =r i /R i and we refer toλ as the target homeostatic stretch. From the Hagen-Poiseuille law, the is available to G&R through WSS signaling from endothelial cells. Hence, the time window average phenotypic modulation stimulus, C, of the vascular SM depends only on this known ratio λ/λ, giving Dividing both sides by λ gives where the effect of the growth history is captured by a SM area-weighted harmonic mean of the stretch: The SM work density can then be formally defined as where the lower limit of the integration is arbitrary owing to the constant term of the primitive function. In a homeostatic state, taking the limit t → ∞ for Equation (17) gives λ m ≡λ, so that the SM work density in the homeostatic state becomeŝ The derivatives of the SM work density (i.e.,Ŵ a ), according to the constitutive relation, give the active force in pathological conditions, which can then be normalized to its homeostatic value obtained from a normal hemodynamics. Here, it was assumed that impaired active response is linearly correlated to the phenotypic shift of SM from the contractile to synthesis state. Hence, impairment was quantified in the G&R model through the ratio of CoA to homeostatic active stress, named the impairment index, and compared to corresponding values obtained empirically using wire myography, i.e., percentage impaired contractile or relaxation behavior of CoA rabbits relative to the control group.
It is worth mentioning that, since the hemodynamic modulation stimulus (C λ λ ) is an unknown function, it was represented using a truncated Taylor expansion around the homeostatic condition, λ λ = 1.
whereĈ and C are nondimensional constants representing the phenotypic modulation stimulus of SM under homeostatic and pathological conditions, respectively. Under homeostasis, S = S(C[1]) = S(Ĉ), which, according to Equation (14), givesĈ = √ − ln 1 −Ŝ/S M = 0.8326 for our choice of parameters. Therefore, we assumedĈ as constant and tuned C in the G&R model to replicate phenotypic modulation and associated active vascular impairment obtained empirically using an impairment index defined in myography results as detailed below.

Strain Energy Density Function for Vessel Constituents
Currently, we define a general form of a strain energy density function to describe the passive action of vessel constituents. In general, if it is assumed only one family of collagen and SM is circumferentially oriented, the strain energy density function for both can be defined similar to elastin using a hybrid neo-Hookean and Fung-type material [30]: where c k i , i ∈ (1, 2, 3) are constitutive parameters characterizing material properties of constituent k ∈ (c, m, e). Now, we have all strain energy and SM work densities in terms of λ, which is a function of current and generation time, i.e., λ(t, τ). Hence, the governing equation of G&R can be derived and discretized for the numerical simulation described in the next section.

Simulation of Stress Mediated G&R
For a pressurized thin-walled incompressible elastic tube with radius r and thickness h, Hoop stress can be derived by force equilibrium in the radial direction: where h r i and denote wall thickness and internal radius, respectively. For laminar flow of a Newtonian fluid with viscosity, η, and velocity, u, the Hagen-Poiseuille equation describes WSS as Considering that, for a thin-walled tube, area can be described as A ∼ = 2πr i h, rewriting Hoop stress and Hagen-Poiseuille in terms of A and dividing by the hemostatic state denoted by hat sign (ˆ), one can derive the current to homeostatic area ratio in terms of mechanical stresses acting on the wall: This equation demonstrates how two known local stimuli, namely, the WSS and IWS, mediate to reproduce the ratioÂ/A and, thus, the appropriate direction of growth: increased rate of production forÂ/A > 1 and decreased rate of production forÂ/A < 1. Taylor series expansion of Equation (24) around the homeostatic stateÂ/A = 1 giveŝ Stress-mediated remodeling theory identifies K σ =σ −1 θ and K φ = 2 3φ −1 , and the basal production rate analogous to the definition of mass production rate as a linear function of mechanical stimulus is as follows [31,32]: with m k 0 being the basal mass production rate of constituent k ∈ (c, m, e). In summary, Equation (26) governs the evolution of the vessel constituents through Equation (1) where the mass per unit area for each vessel constituent is defined. WSS can be derived either from the Hagen-Poiseuille law or from computational fluid dynamics (CFD) using a reduced-order 1D model (ROM). The intramural stress is quantified from constitutive relations and strain energy, and SM work densities are defined according to the growth history in the previous section, 0 ≤ τ ≤ t. Since the right-hand side of Equation (26) depends on the growth history, Equation (1) is the integral form of a delay differential equation (DDE) that is solved using an explicit numerical scheme described below.
The distribution of remaining mass per unit area is defined as such that M k (t, τ)dτ represents the amount of constituent k that remains at time t, and produced within the time interval [τ, τ + dτ). It follows that for any constant time-step ∆t, and that Time is discretized using time-step ∆t, with the following discretization points: Then, M k becomes discretized as From Equation (28), we also have An explicit numerical scheme is used to integrate the evolution of the vessel constituents using Equation (1) with Equation (26). Slowly varying quantities, such as m k (t) and λ m , are treated explicitly, while λ, which changes instantaneously due to changes in the controlling parametersλ(t) and p(t), is treated implicitly. Therefore, integrals are approximated using a trapezoidal rule throughout. For a compact notation, we define a summation symbol used for the trapezoidal rule: The equilibrium Equations (7) and (16) give This equation is solved numerically for λ n+1 withλ n+1 as an initial guess using the MATLAB numerical solver (MathWorks, Natick, MA, USA). If the solution is not unique, the closest solution toλ n+1 is chosen [26].
The stress measure σ n+1 is computed using stress difference terms quantified through differentiating strain energy or SM work density with respect to λ and then multiplying by λ: For each constituent k ∈ (c, m, e), the newly produced materials are computed according to Equations (26) and (29): Moreover, the remaining materials M k n+1,j , j = 0, 1, . . . , n from production at earlier times is computed using Equation (32). The mass per unit area of each constituent is calculated from Equation (1) using the trapezoidal rule: The weighted harmonic mean of the stretch is also expressed through Equation (17) and the trapezoidal rule: Using the trapezoidal rule and Equation (9) for one constituent gives which is a function of λ to be used for implicit calculation of λ n+2 in the subsequent time-step. With these operations, all variables are obtained for the timestep (n + 1)∆t, and the new functions ∂ ∂λ M k W k are defined for use in the following timestep. Here, G k represents the pre-stretch for constituent k and assumed to be constant and equal to 1.08 and 1.4 for collagen and elastin, respectively [28].

G&R Model Fitting
Constitutive parameters of the derived G&R model were tuned to simulate the kinetics of growth and remodeling in response to CoA-induced mechanical stimuli. More specifically, results of the simulations for arterial thickening, stiffening, and dysfunction were compared to measurements of these precursors of hypertension quantified among nine CoA groups (details of empirical measurements are provided in the next section). Constitutive parameters were tuned to minimize the following weighted objective function: where SEE k is the sum of square error of the model vs. empirically measured hypertension precursors with k {sti f f ness, thickness, dys f ucntion}, and w k represents the weight constant for the k-th hypertension precursor. A genetic algorithm was used in MATLAB [33] to solve the nonlinear optimization problem by setting the population size of 1000 and function tolerance of 1 × 10 −7 . Collectively, the resulting fit represented a mathematical model tuned to predict thickness evolution (from longitudinal Doppler ultrasound images) and terminal stiffness and dysfunction (from invasive final material testing and myograph analysis). Table 1 lists the G&R parameters in the model with references and/or values, where applicable.

Empirical Quantification of Hypertension Precursors
Under pathological conditions, G&R can be described through SM synthesis and collagen turnover [16] that are both assumed to be stress-driven [31] and can, therefore, be correlated to the CoA-induced mechanical stimuli. Identifying associated G&R constitutive parameters, thus, requires determination of (1) mechanical stimuli, and (2) arterial dysfunc-tion for an array of cases with different severity and durations of CoA. For this purpose, a rabbit model of CoA accommodating various severities and durations of the mechanical stimuli was implemented [8]. Hypertension precursors including arterial thickening, stiffening, and dysfunction were quantified using an array of noninvasive and invasive experimental protocols [8]. Hypertension precursors included (1) arterial thickening, (2) stiffening, and (3) dysfunction, which were quantified following previously published protocols for a rabbit model mimicking human CoA. Figure 1 summarizes the workflow for the experimental protocol and quantities measured using each modality. After IACUC approval, CoA was surgically induced under isoflurane anesthesia in New Zealand white rabbits (total 63; n = 5-11/group) by tying suture to varying diameters around the aorta at an age of~10 weeks [8]. Tying the suture at different diameters, i.e., 2.06, 1.63, and 1.40 mm, resulted in various clinically important severities named mild (peak-to-peak CoA gradient ≤13 mmHg), intermediate (peak-to-peak CoA gradient 13-20 mmHg), and severe (peak-to-peak CoA gradient ≥20 mmHg). Importantly, the putative threshold suggestive of treatment for CoA in humans is ≥20 mmHg [37].
To investigate the effect of duration of the mechanical stimuli caused by the CoA, sutures with different dissolving properties, i.e., rapid dissolvable, dissolvable, and permanent, were used to present stenoses for different durations of~2, 5, and 22 weeks named short, long, and prolonged CoA, respectively. Importantly, sutures were tied around the proximal descending thoracic aorta (distal to the left subclavian branch) where CoA most often presents. Together, this combination resulted in nine study groups (three severities and three durations) plus a control group.

Temporal Monitoring of Morphology and Hemodynamics
Using body weight at arrival and according to the growth rate for New Zealand White rabbits [38], age was estimated and then scaled to that of humans according to the duration of different phases of life [39] (Figure 1). Temporal hemodynamic and morphological evaluation was performed as a function of age using Doppler ultrasound imaging following protocols similar to human transthoracic echocardiography. Specifically, Doppler color flow and spectral Doppler images ( Figure 2A) were used to identify pre-and post-stenotic flow and associated peak velocity for Doppler-based hemodynamic assessment, while b-mode imaging ( Figure 2B) was used to quantify morphological properties, such as aorta diameter, wall thickness, and percentage area obstruction, using the protocols detailed elsewhere [8,9]. Moreover, phase-contrast magnetic resonance imaging (PC-MRI) was performed to characterize morphological properties, as well as quantify cardiac output and flow distribution to the main aortic branches [8,9]. Lastly, high-fidelity BP measurements were performed at the end of the protocol via invasive catheterization before harvesting of aortic tissue for further invasive characterization, e.g., material properties through uniaxial extension testing and vascular dysfunction through myograph analysis ( Figure 1). Together, these empirical measurements were quantified as surrogates (i.e., precursors) of hypertension including arterial thickening, stiffening, and dysfunction, as detailed in the next subsections. . Experimental protocol to characterize hypertension precursors in response to CoA in rabbits. Male New Zealand White rabbits were randomly selected for a control group or designated to undergo descending thoracic CoA (nine study groups: three severities and three durations) via left thoracotomy in the third intercostal space. Weekly Doppler ultrasound imaging characterized local hemodynamic and morphology adaptation to CoA before PC-MRI prior to the final week. At the end of the protocol, invasive characterization was conducted for blood pressure measurement, and the aorta was dissected to quantify arterial stiffening and dysfunction through uniaxial extension testing and myograph analysis, respectively.

Temporal Monitoring of Morphology and Hemodynamics
Using body weight at arrival and according to the growth rate for New Zealand White rabbits [38], age was estimated and then scaled to that of humans according to the duration of different phases of life [39] (Figure 1). Temporal hemodynamic and morphological evaluation was performed as a function of age using Doppler ultrasound imaging following protocols similar to human transthoracic echocardiography. Specifically, Doppler color flow and spectral Doppler images ( Figure 2A) were used to identify pre-and post-stenotic flow and associated peak velocity for Doppler-based hemodynamic assessment, while b-mode imaging ( Figure 2B) was used to quantify morphological properties, such as aorta diameter, wall thickness, and percentage area obstruction, using the protocols detailed elsewhere [8,9]. Moreover, phase-contrast magnetic resonance imaging (PC-MRI) was performed to characterize morphological properties, as well as quantify cardiac output and flow distribution to the main aortic branches [8,9]. Lastly, high-fidelity BP measurements were performed at the end of the protocol via invasive catheterization before harvesting of aortic tissue for further invasive characterization, e.g., material properties through uniaxial extension testing and vascular dysfunction through myograph analysis ( Figure 1). Together, these empirical measurements were quantified as surrogates (i.e., precursors) of hypertension including arterial thickening, stiffening, and dysfunction, as detailed in the next subsections. Experimental protocol to characterize hypertension precursors in response to CoA in rabbits. Male New Zealand White rabbits were randomly selected for a control group or designated to undergo descending thoracic CoA (nine study groups: three severities and three durations) via left thoracotomy in the third intercostal space. Weekly Doppler ultrasound imaging characterized local hemodynamic and morphology adaptation to CoA before PC-MRI prior to the final week. At the end of the protocol, invasive characterization was conducted for blood pressure measurement, and the aorta was dissected to quantify arterial stiffening and dysfunction through uniaxial extension testing and myograph analysis, respectively.

Vascular Stiffening
To investigate coarctation-induced arterial stiffening, material properties were quantified in the proximal descending thoracic aorta (region between the left subclavian artery and the coarctation). This region was selected due to its exposure to adversely

Vascular Stiffening
To investigate coarctation-induced arterial stiffening, material properties were quantified in the proximal descending thoracic aorta (region between the left subclavian artery and the coarctation). This region was selected due to its exposure to adversely elevated coarctation-induced BP. Material properties were characterized through uniaxial extension testing (MTS Criterion Load Frame, MTS, Minneapolis, MN, USA) at 37 • C in an environmental chamber (MTS Bionix EnviroBath, Minneapolis, MN, USA). Tissues were dissected in the circumferential direction with a length-to-width ratio of~2.6 and preconditioned by stretching to 10% of the gauge length. Extension testing was performed at 10 mm/min until hyperelastic behavior appeared. Results were interpreted in stressstrain curves to be compared to G&R model predictions quantified via Equation (34) at the experimental endpoint.

Vascular Thickening
Doppler b-mode images were used to quantify aortic wall thickness in the proximal descending thoracic aorta ( Figure 2B). All thickness measurements were made in triplicate with the mean value reported for each date. Body weight was also measured weekly using a scientific scale. Lastly, the thickness evolution was normalized to interpolated body weights at each timepoint assuming scaling of the body weight according to a sigmoid function ( Figure 3). The sigmoid function was selected specifically to allow for a relatively low initial growth rate during the gestation period (~32 days) assuming an average of 50 g for newborn rabbit weight [39], as well as stable body weight, upon reaching adulthood. Doppler ultrasound follow-up was performed weekly to quantify the temporal evolution of hemodynamic changes (A, spectral Doppler) and morphology (B, Doppler b-mode) in response to the CoA.

Vascular Stiffening
To investigate coarctation-induced arterial stiffening, material properties were quantified in the proximal descending thoracic aorta (region between the left subclavian artery and the coarctation). This region was selected due to its exposure to adversely elevated coarctation-induced BP. Material properties were characterized through uniaxial extension testing (MTS Criterion Load Frame, MTS, Minneapolis, MN, USA) at 37 °C in an environmental chamber (MTS Bionix EnviroBath, Minneapolis, MN, USA). Tissues were dissected in the circumferential direction with a length-to-width ratio of ~2.6 and preconditioned by stretching to 10% of the gauge length. Extension testing was performed at 10 mm/min until hyperelastic behavior appeared. Results were interpreted in stressstrain curves to be compared to G&R model predictions quantified via Equation (34) at the experimental endpoint.

Vascular Thickening
Doppler b-mode images were used to quantify aortic wall thickness in the proximal descending thoracic aorta ( Figure 2B). All thickness measurements were made in triplicate with the mean value reported for each date. Body weight was also measured weekly using a scientific scale. Lastly, the thickness evolution was normalized to interpolated body weights at each timepoint assuming scaling of the body weight according to a sigmoid function ( Figure 3). The sigmoid function was selected specifically to allow for a relatively low initial growth rate during the gestation period (~32 days) assuming an average of 50 g for newborn rabbit weight [39], as well as stable body weight, upon reaching adulthood.  . Representative thickness and body weight evolution over time for a CoA rabbit. Thickness (squares) was quantified as a function of age and through weekly ultrasound imaging at follow-up where 2D sagittal sections of the aorta were identified in b-mode images, and triplicate thickness measurements were performed with mean values shown at each timepoint. Body weight evolution was extrapolated using a sigmoid fit (dashed line) to the measured body weights (circles).

Vascular Dysfunction
Vascular specimens (3-4 mm rings) were carefully sectioned from the proximal descending thoracic aorta and assessed for vascular dysfunction using a protocol described previously [8]. Briefly, experimental assessment of active vascular dysfunction was conducted to observe SM contraction via phenylephrine (PE) in a half-log increasing dose responses from 10 −9 to 10 −5 . To assess endothelium-dependent relaxation via nitric oxide (NO), arteries were precontracted with PE to the EC 50 concentration. Once plateau was achieved, cumulative addition of the endothelium-dependent muscarinic receptor agonist acetylcholine (ACh; 10 −9 to 10 −5 M) was initiated, and relaxation response curves were quantified as a percentage of precontracted active tension. Area under the dose response curves (AUC) were then quantified as an aggregate measure of vasoactive dysfunction and normalized to that of the control group in a normalized quantity we refer to as the impairment index. All quantifications were performed in duplicate for paired channels with mean values reported for each rabbit. The impairment index was used as a quantitative nondimensional measure of active vascular impairment due to shift in SM phenotype to the dedifferentiated state and underlying endothelial dysfunction. This index was also quantified computationally using the CoA-to-control SM active stress ratio (Equation (18)) and used as a tool to tune the phenotypic modulation stimuli constant C in the G&R model.
2.5.6. WSS and IWS Evolution WSS (φ) and IWS (σ) are among the most important stress components deriving the rate of stress-mediated growth and remodeling as defined by Equation (26). WSS evolution was calculated on the basis of empirical data from the proximal descending thoracic aorta by implementing peak Doppler velocity readings in the Hagen-Poiseuille law assuming a parabolic velocity profile (Equation (23)). Additionally, IWS in the circumferential direction was quantified from the force equilibrium (Equation (22)). The radius and thickness in this equation were quantified from Doppler b-mode images as described earlier. As mentioned above, mean arterial pressure, i.e., MAP (p), was quantified through high-fidelity BP measurements using catheterization and described as a linear function of percentage area obstruction ( Figure 4) determined from Doppler ultrasound images. oxide (NO), arteries were precontracted with PE to the EC50 concentration. Once plateau was achieved, cumulative addition of the endothelium-dependent muscarinic receptor agonist acetylcholine (ACh; 10 −9 to 10 −5 M) was initiated, and relaxation response curves were quantified as a percentage of precontracted active tension. Area under the dose response curves (AUC) were then quantified as an aggregate measure of vasoactive dysfunction and normalized to that of the control group in a normalized quantity we refer to as the impairment index. All quantifications were performed in duplicate for paired channels with mean values reported for each rabbit. The impairment index was used as a quantitative nondimensional measure of active vascular impairment due to shift in SM phenotype to the dedifferentiated state and underlying endothelial dysfunction. This index was also quantified computationally using the CoA-to-control SM active stress ratio (equation 18) and used as a tool to tune the phenotypic modulation stimuli constant ′ in the G&R model.

WSS and IWS Evolution
WSS ( ) and IWS ( ) are among the most important stress components deriving the rate of stress-mediated growth and remodeling as defined by Equation (26). WSS evolution was calculated on the basis of empirical data from the proximal descending thoracic aorta by implementing peak Doppler velocity readings in the Hagen-Poiseuille law assuming a parabolic velocity profile (Equation (23)). Additionally, IWS in the circumferential direction was quantified from the force equilibrium (Equation (22)). The radius and thickness in this equation were quantified from Doppler b-mode images as described earlier. As mentioned above, mean arterial pressure, i.e., MAP ( ), was quantified through high-fidelity BP measurements using catheterization and described as a linear function of percentage area obstruction ( Figure 4) determined from Doppler ultrasound images.

Statistical Analysis
Descriptive statistics are presented for continuous variables as the mean ± standard error of the mean (SEM). Unbalanced one-way analysis of variance (ANOVA) was used to assess significant differences across groups using Dunnett's multiple comparison post hoc analysis with 5% level of significance in GraphPad Prism version 9.5.1, GraphPad Software, San Diego, CA, USA, www.graphpad.com (accessed on 12 October 2022). Pearson's correlation and linear regression analysis examined the relationships and r-squared for goodness-of-fit. G&R model prediction was validated through impairment index values obtained from rabbits not used for fitting and in terms of percentage error and R-squared. D'Agostino and Pearson, Anderson-Darling, Shapiro-Wilk, and Kolmogorov-Smirnov tests were all used to assess the normality of prediction errors at the 0.05 significance level.

Results
Following the use of permanent suture in the prolonged CoA groups, the stenosis was present for~22 weeks, resulting in adversely elevated BP. In contrast, degradation of the suture in the short and long CoA groups restored aortic diameter and, consequently, also restored associated elevation in BP. Area obstruction levels vary as a function of experimental group as shown in Figure 5 due to rabbit growth and resorption of dissolvable sutures.
Pearson's correlation and linear regression analysis examined the relationships and rsquared for goodness-of-fit. G&R model prediction was validated through impairment index values obtained from rabbits not used for fitting and in terms of percentage error and R-squared. D'Agostino and Pearson, Anderson-Darling, Shapiro-Wilk, and Kolmogorov-Smirnov tests were all used to assess the normality of prediction errors at the 0.05 significance level.

Results
Following the use of permanent suture in the prolonged CoA groups, the stenosis was present for ~22 weeks, resulting in adversely elevated BP. In contrast, degradation of the suture in the short and long CoA groups restored aortic diameter and, consequently, also restored associated elevation in BP. Area obstruction levels vary as a function of experimental group as shown in Figure 5 due to rabbit growth and resorption of dissolvable sutures. Figure 5. Percentage area obstruction vs. days after CoA surgery. CoA morphology was evaluated on a weekly basis noninvasively using ultrasound imaging. Proximal aorta and CoA diameters were quantified in triplicate, and the average value was plotted (scatter) at each date. The percentage area obstruction (%AO = 100 × (Dprox 2 − DCoA 2 )/Dprox 2 ) was, therefore, quantified and plotted vs. days after surgery. For the dissolving sutures (i.e., short and long duration CoA groups), a plateau followed by exponential decay function was fitted to the %AO, while, for the permanent suture (i.e., prolonged CoA group), an exponential function was used to fit the data. days after surgery. For the dissolving sutures (i.e., short and long duration CoA groups), a plateau followed by exponential decay function was fitted to the %AO, while, for the permanent suture (i.e., prolonged CoA group), an exponential function was used to fit the data.
Arterial thickening. Figure 6 shows the normalized thickness evolution for the CoA vs. control groups over the course of the protocol. Normalized thickness initially decreased in all rabbit groups and eventually plateaued with body weight after~150 days of age. Overall, arterial thickening was observed in intermediate to severe CoA (peak-to-peak CoA gradient ≥13 mmHg) when presented for long (~5 weeks) or prolonged durations (~22 weeks). The pattern of arterial thickening was similar to the control group for short presentation of the CoA (i.e., all severities for~2 weeks). In the mild CoA group (peakto-peak CoA gradient <13 mmHg); however, thickness for the prolonged CoA group was slightly larger than that of the control at the plateau. After parameter calibration, the G&R model successfully replicated thickness evolutions. Table 2 reports the parameters of the G&R model obtained after the optimization procedure. The adjusted R-squared values for the resulting thickness evolution fits were in the range of 0.98 to 1.00 for all groups. presentation of the CoA (i.e., all severities for ~2 weeks). In the mild CoA group (peak-topeak CoA gradient <13 mmHg); however, thickness for the prolonged CoA group was slightly larger than that of the control at the plateau. After parameter calibration, the G&R model successfully replicated thickness evolutions. Table 2 reports the parameters of the G&R model obtained after the optimization procedure. The adjusted R-squared values for the resulting thickness evolution fits were in the range of 0.98 to 1.00 for all groups. Coarctation-induced arterial thickening (shades of gray) was quantified from Doppler ultrasound b-mode images. To compare with control group, thicknesses were normalized to body weight (black). CoA groups included all combinations of severity (mild, intermediate, and severe) and duration (short, long, and prolonged) studied. G&R model constitutive parameters were tuned to fit empirical measurements (n = 2/group) with adjusted R 2 characterizing goodness of fit.  Figure 6. Thickness evolution for CoA rabbits (scatter) and G&R model predictions (lines). Coarctation-induced arterial thickening (shades of gray) was quantified from Doppler ultrasound b-mode images. To compare with control group, thicknesses were normalized to body weight (black). CoA groups included all combinations of severity (mild, intermediate, and severe) and duration (short, long, and prolonged) studied. G&R model constitutive parameters were tuned to fit empirical measurements (n = 2/group) with adjusted R 2 characterizing goodness of fit. * Collagen and smooth muscle were described by an exponential strain energy density function [30] that requires c 1 = 0. ** Elastin was described by a neo-Hookean [17] strain energy density function that requires c 2 = 0 and also drops c 3 from Equation (21). *** c 1 and c 2 are in kPa, and c 3 is nondimensional. Values represent the average ± SD. † Pressures are in units of mmHg.
Arterial stiffening. Figure 7 represents the empirical stress-strain curves quantified in the proximal descending thoracic aorta exposed to adversely elevated BP in CoA groups. Results from each CoA group are compared to the control group in each plot. Overall, the pattern of arterial stiffening was observed by increased stress at the same stretch ratio for groups exposed to a prolonged CoA (~22 weeks). This was also observed in the long CoA groups (~5 weeks), but only for the intermediate and severe CoA (i.e., peak-to-peak CoA gradient ≥13 mmHg). To investigate the rate of stress-mediated remodeling and its effect on the vessel stiffness, constitutive parameters for collagen, SM, and elastin were tuned to fit the empirically measured stress-stretch curves at the end of the protocol (Figure 7, red lines). The constitutive material parameters obtained for c 1 , c 2 , and c 3 are listed in Table 2.
Arterial stiffening. Figure 7 represents the empirical stress-strain curves quantified in the proximal descending thoracic aorta exposed to adversely elevated BP in CoA groups. Results from each CoA group are compared to the control group in each plot. Overall, the pattern of arterial stiffening was observed by increased stress at the same stretch ratio for groups exposed to a prolonged CoA (~22 weeks). This was also observed in the long CoA groups (~5 weeks), but only for the intermediate and severe CoA (i.e., peak-to-peak CoA gradient ≥13 mmHg). To investigate the rate of stress-mediated remodeling and its effect on the vessel stiffness, constitutive parameters for collagen, SM, and elastin were tuned to fit the empirically measured stress-stretch curves at the end of the protocol (Figure 7, red lines). The constitutive material parameters obtained for c1, c2, and c3 are listed in Table 2. Figure 7. Empirical stress-stretch curves for CoA (gray) vs. control (black) groups and G&R model prediction (red). Error bars representing the SEM were quantified and averaged at each strain energy level discretized over the range of 0 to 3 MPa. The constitutive material parameters identified are listed in Table 2.
Vascular dysfunction. Overall, CoA rabbits showed impaired vascular response to halflog increasing dose of PE and ACh as compared to the control group. The extent of vascular impairment showed patterns correlating with both severity and duration of the CoA-induced mechanical stimuli. Figure 8 shows the normalized effective force response and associated AUC ( Figure 8D,H,L-O) for CoA vs. control groups. The normalized effective force showed an intact logistic contraction response in the control rabbits peaking at 1.24 ± 0.07, whereas CoA groups peaked to smaller normalized tension values ranging from 0.90 ± 0.06 to 1.09 ± 0.13. This diminished contractile capacity was also observed through AUC among CoA rabbits ( Figure 8D,H,L-O). For example, when CoA was at the most severe level, i.e., peak-to-peak CoA gradient ≥20 mmHg, the AUC was significantly Figure 7. Empirical stress-stretch curves for CoA (gray) vs. control (black) groups and G&R model prediction (red). Error bars representing the SEM were quantified and averaged at each strain energy level discretized over the range of 0 to 3 MPa. The constitutive material parameters identified are listed in Table 2.
Vascular dysfunction. Overall, CoA rabbits showed impaired vascular response to halflog increasing dose of PE and ACh as compared to the control group. The extent of vascular impairment showed patterns correlating with both severity and duration of the CoA-induced mechanical stimuli. Figure 8 shows the normalized effective force response and associated AUC ( Figure 8D,H,L-O) for CoA vs. control groups. The normalized effective force showed an intact logistic contraction response in the control rabbits peaking at 1.24 ± 0.07, whereas CoA groups peaked to smaller normalized tension values ranging from 0.90 ± 0.06 to 1.09 ± 0.13. This diminished contractile capacity was also observed through AUC among CoA rabbits ( Figure 8D,H,L-O). For example, when CoA was at the most severe level, i.e., peak-to-peak CoA gradient ≥20 mmHg, the AUC was significantly smaller than in control rabbits, regardless of duration ( Figure 8L). Similarly, when CoA was present for the prolonged period (~22 weeks), significant impairment was observed in AUC, regardless of severity ( Figure 8O). A nonsignificant trend toward impaired PE contraction was observed in all other groups. Figure 9 shows arterial relaxation curves in response to the endothelium-dependent agonist ACh for the descending thoracic aorta proximal to the CoA (exposed to adversely elevated coarctation-induced BP). Overall, control rabbits (black, n = 11) showed intact endothelial function with relaxation to 53% ± 4% of the precontraction force at 10 −7 M ACh concentration. However, this vasodilatory action of ACh was impaired in CoA rabbits (gray, n = 5-10). The extent of impaired ACh relaxation correlated with the severity of CoA-induced mechanical stimuli. More specifically, significant decreases in AUC ( Figure 9D,H,L-O) were observed in the severe CoA group (peak-to-peak CoA gradient ≥20 mmHg) regardless of the duration of the CoA ( Figure 9L). Consequently, the quantified impairment index (CoA-to-control AUC ratio) was 0.52 ± 0.27, 0.55 ± 0.13, and 0.52 ± 0.16 for short, long, and prolonged duration groups, respectively. Rabbits with intermediate CoAs also developed impaired ACh relaxation response that was significant when exposed for >2 weeks, i.e., long and prolonged duration groups ( Figure 9N). Associated impairment indices were 0.67 ± 0.19 and 0.85 ± 0.08, respectively. A trend toward impairment was observed in the mild CoA groups, but differences did not reach significance ( Figure 9M).
Biomedicines 2023, 11,1817 16 of 23 smaller than in control rabbits, regardless of duration ( Figure 8L). Similarly, when CoA was present for the prolonged period (~22 weeks), significant impairment was observed in AUC, regardless of severity ( Figure 8O). A nonsignificant trend toward impaired PE contraction was observed in all other groups. . Error bars represent the SEM. * Significant difference compared to the control group in a two-tailed Student's t-test with a 5% significance level, assuming similar variation among the groups. Figure 9 shows arterial relaxation curves in response to the endothelium-dependent agonist ACh for the descending thoracic aorta proximal to the CoA (exposed to adversely elevated coarctation-induced BP). Overall, control rabbits (black, n = 11) showed intact endothelial function with relaxation to 53% ± 4% of the precontraction force at 10 −7 M ACh concentration. However, this vasodilatory action of ACh was impaired in CoA rabbits (gray, n = 5-10). The extent of impaired ACh relaxation correlated with the severity of CoA-induced mechanical stimuli. More specifically, significant decreases in AUC ( Figure  9D,H,L-O) were observed in the severe CoA group (peak-to-peak CoA gradient ≥20 mmHg) regardless of the duration of the CoA ( Figure 9L). Consequently, the quantified impairment index (CoA-to-control AUC ratio) was 0.52 ± 0.27, 0.55 ± 0.13, and 0.52 ± 0.16 for short, long, and prolonged duration groups, respectively. Rabbits with intermediate CoAs also developed impaired ACh relaxation response that was significant when exposed for >2 weeks, i.e., long and prolonged duration groups ( Figure 9N). Associated impairment indices were 0.67 ± 0.19 and 0.85 ± 0.08, respectively. A trend toward impairment was observed in the mild CoA groups, but differences did not reach significance ( Figure 9M). Phenotypic modulation was characterized computationally through its stimulus term (C λ /λ ) for SM as a function of the current to homeostasis (i.e., control) stretch ratio ( λ /λ ). As indicated in Equation (20),Ĉ and C are the two nondimensional constants that determine the phenotypic modulation in SM proliferation and impaired active response. Figure 10 shows active stress per unit relaxed SM normalized to that of the control group (impairment index) and compared to the corresponding empirical measurements of PE and ACh myograph results. Overall, PE contraction and ACh relaxation were both impaired with indices ranging from 0.55 to 0.94 (PE), and from 0.52 to 0.93 (Ach), as shown in Figure 10. Interestingly, this index correlated with the severity and duration of the mechanical stimuli caused by the CoA with values as low as 0.52 observed in the most severe and prolonged CoA groups. Phenotypic modulation stimulus constitutive parameters (C ) tuned in the G&R model (n = 2) to replicate impairment quantified from wire myography are summarized for each group in Table 3. Phenotypic modulation was characterized computationally through its stimulus term ( [ / ]) for SM as a function of the current to homeostasis (i.e., control) stretch ratio ([ / ]). As indicated in Equation (20), and are the two nondimensional constants that determine the phenotypic modulation in SM proliferation and impaired active response. Figure 10 shows active stress per unit relaxed SM normalized to that of the control group (impairment index) and compared to the corresponding empirical measurements of PE and ACh myograph results. Overall, PE contraction and ACh relaxation were both impaired with indices ranging from 0.55 to 0.94 (PE), and from 0.52 to 0.93 (Ach), as shown in Figure 10. Interestingly, this index correlated with the severity and duration of the mechanical stimuli caused by the CoA with values as low as 0.52 observed in the most severe and prolonged CoA groups. Phenotypic modulation stimulus constitutive parameters ( ) tuned in the G&R model (n = 2) to replicate impairment quantified from wire myography are summarized for each group in Table 3. Dose-response curves were quantified via wire myography and a half-log increase in phenylephrine (PE, gray) and acetylcholine (Ach, pink) to characterize NO active response in SM contraction and relaxation, respectively. A similar index was also quantified by taking the ratio of active stress per unit relaxed SM for CoA relative to control groups using the G&R framework (green) for two rabbits per group. G&R constitutive parameters were tuned to achieve minimum error to empirically quantified active response, i.e., PE and Ach. No significant differences were observed between results from empirical and G&R through two-tailed Student's t-test with a 5% significance level.  Dose-response curves were quantified via wire myography and a half-log increase in phenylephrine (PE, gray) and acetylcholine (Ach, pink) to characterize NO active response in SM contraction and relaxation, respectively. A similar index was also quantified by taking the ratio of active stress per unit relaxed SM for CoA relative to control groups using the G&R framework (green) for two rabbits per group. G&R constitutive parameters were tuned to achieve minimum error to empirically quantified active response, i.e., PE and Ach. No significant differences were observed between results from empirical and G&R through two-tailed Student's t-test with a 5% significance level. For example, Figure 11A shows the impairment index in the clinically important range of severity, as well as critical age phases. The lines of constant impairment index shown can be used as a map for CoA assessment considering both severity of the narrowing and rabbit age. Particularly, an active response <50%, i.e., impairment index <0.50, indicates nonviable tissue as reported by Lopez et al. [45]. Using this criterion, our results indicate that a peakto-peak CoA gradient of >35 mmHg, even when exposed shortly, can permanently diminish the active arterial response ( Figure 11A). The current model predictions show that such functional impairment is expected for a variety of CoA severity and durations including a short exposure to >30 mmHg peak-to-peak CoA gradient, as well as prolonged exposure to a milder CoA gradient (~10 mmHg). Interestingly, in ages younger than the rabbits studied, i.e., log 10 (days) < 1.75, i.e.,~70 rabbit days of age (or~1 human year), model extrapolation showed arterial impairment being strongly correlated to severity of the CoA gradient, thus emphasizing the importance of early intervention to prevent impaired vascular function. This finding further emphasizes the application of such predictive models in clinical decision making, where both severity and duration of the coarctation-induced mechanical stimuli contribute to hypertension precursors developing in the aorta and major branches. CoA severity, shown on the vertical axis, is characterized through peak-to-peak trans-coarctation BP gradient (BPG), whereas duration is represented on a logarithmic scale of age in days. (B) Model error relative to measured impairment index, i.e., area under the active dose-response curves (AUC) for each CoA group normalized to that of the control. Dose-response curves were quantified via wire myography and half-log increase of phenylephrine (PE, gray) and acetylcholine (Ach, pink) to characterize NO active response in SM contraction and relaxation, respectively. Model error was then quantified as a percentage relative to the empirical PE and Ach impairment indices. Rabbits used in the tuning process were excluded for model prediction error.
The current results should be interpreted relative to several potential limitations. In particular, high-fidelity hemodynamic simulation is challenging to perform longitudinally as 3D reconstruction of the aorta over the course of the disease requires multiple MRI acquisitions and often computationally heavy [46,47]. We, therefore, quantified WSS and IWS from the Hagen-Poiseuille law and force equilibrium using empirical measurements of local hemodynamics over the course of disease development. Additionally, the large number of unknowns in the formulation derived for G&R requires a large sample size, and that is challenging to achieve due to the time and cost of experimental procedures needed for longitudinal monitoring of hypertension precursors using our preclinical model. Therefore, many of the constitutive parameters were assumed on the basis of the literature followed by tuning of the stress-mediating G&R parameters (Table 3) to replicate empirical results. It is also worth mentioning that the pathophysiological response of the aorta is also affected by concomitant morphological anomalies including the bicuspid aortic valve, aortic arch/isthmus hypoplasia, and longsegment coarctation. The associated effect on local hemodynamic and aortic remodeling remains to be investigated in future work. impairment index quantified by taking the ratio of active stress per unit relaxed SM for CoA relative to control using the G&R framework. Contour lines indicate regions of constant impairment indices. CoA severity, shown on the vertical axis, is characterized through peak-to-peak trans-coarctation BP gradient (BPG), whereas duration is represented on a logarithmic scale of age in days. (B) Model error relative to measured impairment index, i.e., area under the active dose-response curves (AUC) for each CoA group normalized to that of the control. Dose-response curves were quantified via wire myography and half-log increase of phenylephrine (PE, gray) and acetylcholine (Ach, pink) to characterize NO active response in SM contraction and relaxation, respectively. Model error was then quantified as a percentage relative to the empirical PE and Ach impairment indices. Rabbits used in the tuning process were excluded for model prediction error. Figure 11B shows active response prediction errors for impairment indices representing NO active response in SM contraction and relaxation. Variability in prediction error was observed which was mainly considered to be inherited from biological variability of the impairment index reported in Figure 10. Statistical assessment of errors showed a normal distribution of across all study groups for both contraction and relaxation impairment indices using four normality tests mentioned in the methods section with p-values > 0.1.
The current results should be interpreted relative to several potential limitations. In particular, high-fidelity hemodynamic simulation is challenging to perform longitudinally as 3D reconstruction of the aorta over the course of the disease requires multiple MRI acquisitions and often computationally heavy [46,47]. We, therefore, quantified WSS and IWS from the Hagen-Poiseuille law and force equilibrium using empirical measurements of local hemodynamics over the course of disease development. Additionally, the large number of unknowns in the formulation derived for G&R requires a large sample size, and that is challenging to achieve due to the time and cost of experimental procedures needed for longitudinal monitoring of hypertension precursors using our preclinical model. Therefore, many of the constitutive parameters were assumed on the basis of the literature followed by tuning of the stress-mediating G&R parameters (Table 3) to replicate empirical results. It is also worth mentioning that the pathophysiological response of the aorta is also affected by concomitant morphological anomalies including the bicuspid aortic valve, aortic arch/isthmus hypoplasia, and long-segment coarctation. The associated effect on local hemodynamic and aortic remodeling remains to be investigated in future work.

Conclusions
Results showed that the implemented G&R model could accurately represent hypertension precursors among the studied groups of CoA rabbits. Therefore, the G&R model allows for computational identification of the adverse vascular response as a function of the severity and duration of coarctation-induced stimuli. Clinical translation of this predictive computational model can help in modifying current intervention thresholds and preventing irreversible hypertension precursors that seem to persist after CoA treatment using current clinical guidance.

Conflicts of Interest:
The terms of graduate student funding for Arash Ghorbannia, the corresponding author of the current work, have been reviewed and approved by the Medical College of Wisconsin and Marquette University in accordance with their policies on objectivity in research. The authors have no other financial disclosures or conflict of interest to report.