Application and reduction of a nonlinear hyperelastic wall model capturing ex vivo relationships between fluid pressure, area, and wall thickness in normal and hypertensive murine left pulmonary arteries

Pulmonary hypertension is a cardiovascular disorder manifested by elevated mean arterial blood pressure (>20 mmHg) together with vessel wall stiffening and thickening due to alterations in collagen, elastin, and smooth muscle cells. Hypoxia‐induced (type 3) pulmonary hypertension can be studied in animals exposed to a low oxygen environment for prolonged time periods leading to biomechanical alterations in vessel wall structure. This study introduces a novel approach to formulating a reduced order nonlinear elastic structural wall model for a large pulmonary artery. The model relating blood pressure and area is calibrated using ex vivo measurements of vessel diameter and wall thickness changes, under controlled pressure conditions, in left pulmonary arteries isolated from control and hypertensive mice. A two‐layer, hyperelastic, and anisotropic model incorporating residual stresses is formulated using the Holzapfel–Gasser–Ogden model. Complex relations predicting vessel area and wall thickness with increasing blood pressure are derived and calibrated using the data. Sensitivity analysis, parameter estimation, subset selection, and physical plausibility arguments are used to systematically reduce the 16‐parameter model to one in which a much smaller subset of identifiable parameters is estimated via solution of an inverse problem. Our final reduced one layer model includes a single set of three elastic moduli. Estimated ranges of these parameters demonstrate that nonlinear stiffening is dominated by elastin in the control animals and by collagen in the hypertensive animals. The pressure–area relation developed in this novel manner has potential impact on one‐dimensional fluids network models of vessel wall remodeling in the presence of cardiovascular disease.


Introduction
Pulmonary hypertension (PH), encompassing several cardiovascular disorders and manifested by a mean pulmonary arterial blood pressure (BP) above 20 mmHg, is commonly classified into five disease groups [14,31].One of these, group 3:"pulmonary hypertension due to lung disease" includes patients with PH induced by hypoxia (HPH).This disease type can be studied in mice with PH induced by placing the animals in a low oxygen (hypoxic) environment.The response of the cardiovascular system is stiffening and thickening of the pulmonary arteries accompanied by an increase in BP to PH levels.The aim of this study is to use mathematical modeling to devise a relationship between BP and vessel lumen area that characterizes the structural remodeling of the underlying tissues.For this PH group, vascular remodeling typically starts in the small arteries, proceeding to the large arteries as the disease advances [16,27].The arterial wall comprises three layers, the intima, a single layer of endothelial cells, the media which contains large amounts of elastin and smooth muscle cells, and the adventitia mainly composed of collagen (Fig. 1a).In animal models of group 3 PH, vessels stiffen largely due to collagen accumulation [20,35] and smooth muscle cell proliferation is known to increase the thickness of the vessel wall [36].
One advantage of characterizing how PH impacts the pressure-area relationship is that the resulting model can be incorporated into one-dimensional (1D) fluid dynamics network models used extensively to study hemodynamics in both systemic [1,2,5,12,18,34] and pulmonary [6,21,25] arteries.1D fluid dynamics models are especially well suited to predict flow distribution along the network and wave-propagation, but accurate predictions require appropriate specification of the pressure-area interaction.Moreover, 1D models can be readily calibrated to in vivo geometry, flow and/or BP measurements.The 1D fluid dynamics models are derived from the Navier-Stokes equations combined with a state equation relating blood pressure and vessel area, often formulated using an empirical or simple elastic wall model.These simpler models have the advantage of being specified using a small number of parameters [11,18,24,33], but the disadvantage that the manner in which tissue remodeling associated with disease translates to the model is unclear.While complex tissue mechanics models exist [8,26,38], they have not been integrated with 1D fluid dynamics models.One state-of-the-art tissue mechanics model is the two-layer nonlinear hyperelastic model developed by Holzapfel, Gasser, and Ogden [8] (HGO model) that captures ex vivo biomechanical deformation of the vessel wall.While this model is complex, it includes parameters that more directly and realistically represent structural elements and constituents within the underlying biological soft tissues.
In this study, we formulate and systematically reduce a nonlinear hyperelastic structural wall model for the large pulmonary arteries, generating a novel pressurearea relation that can characterize remodeling in HPH.The model will be calibrated to ex vivo biomechanical deformation and wall thickness measurements from control and hypertensive mice.To do so, we start by formulating a two-layer, anisotropic vessel wall model using the HGO model formulation [8].In addition to anisotropy and multiple layers, this model accounts for residual stresses, known to be significant in large pulmonary arteries as evidenced by a large opening angle arising when rings from excised vessels are cut.The rings are obtained from cutting "a slice" normal to the axial direction, and the opening angle is determined from a radial cut through the ring's circumference [9,37].Based on this approach, complex relations determining the dependence of vessel area and wall thickness on blood pressure (BP) are derived.Our model will be calibrated to ex vivo measurements of vessel diameter and wall thickness as a function of pressure in the left pulmonary artery (LPA) in control (CTL) and hypertensive (HPH) mice [32].Since our model is complex, containing 16 parameters, calibration using data is challenging.To remedy this problem we use sensitivity analysis and subset selection [3,4,15,23] to identify the simplest model and a set of sensitive and identifiable parameters that can be estimated using the model and available data.

Models and Methods
In 1D cardiovascular fluid dynamics network models, the dependent variables are the transmural blood pressure p(z, t) (mmHg) (the difference between blood pressure in the vessel and the surrounding tissue), the vessel lumen area a(z, t) (cm 2 ), and the average flow q(z, t) (cm 3 /s) through the vessel.t (s) denotes time and z (cm) is the axial coordinate (along the length of the vessel).The flow is approximated by integrating over the vessel's cross-section, i.e. q(z, t) = 2π r(z,t) 0 u z r dr, where a(z, t) = π(r(z, t)) 2 , r(z, t) (cm) is the vessel radius, and u z (cm/s) is the axial component of the fluid velocity.This equation is evaluated by specifying the velocity profile across the lumen [34,18].
The 1D fluid dynamics model is derived from an approximation to the Navier-Stokes equations that conserves volume and momentum.The derivation is achieved by assuming that the flow is incompressible and Newtonian, that the vessels are cylindrical, and that the blood is incompressible, viscous, and homogeneous with constant density ρ (g/cm 3 ) and blood viscosity µ (g/(cm s)).Under these assumptions, the 1D fluid dynamics model can be formulated as where ν = µ/ρ (cm 2 /s) is the Newtonian fluid kinematic viscosity and δ (cm) is a boundary layer thickness parameter introduced via the velocity profile equation to ensure no-slip at the vessel walls.The system of equations is closed by introducing a wall model relating the transmural blood pressure p(z, t) and lumen area a(z, t).This study derives such a relation by treating the wall as a hyperelastic material integrating the two layer model by Holzapfel, Gasser and Ogden, often referred to as the HGO model [8].This model incorporates nonlinear effects of residual stresses, anisotropy, material and geometric nonlinearities, and contributions of key wall constituents (collagen and elastin) within the vessel wall layers.A schematic of the wall constituents is shown in Fig. 1, and Table 4 (in the Appendix) lists the model parameters and their units.

Deformation
The model is formulated in terms of three configurations of the vessel wall: (i) a stress-free reference state Ω 0 (Fig. 1b) represented by a continuous arc of a cylindrical ring free of all residual stresses; (ii) an intermediate load-free configuration (not shown) represented by a closed cylindrical ring in the absence of fluid flow; and (iii) a current configuration Ω (Fig. 1c) representing the deformed vessel as fluid flows through the vessel lumen in an ex vivo or in vivo setting.
Specifically, Ω 0 approximates the process of excising a vessel segment, extracting a cross-section approximated as a thin cylindrical ring, and then making a single radial cut along the ring's circumference.It is denoted by where (R, Θ, Z) are Lagragian cylindrical (polar) coordinates, α is the opening angle, L is the reference axial length, and R in and R out are the inner and outer radii, respectively.Similarly, the current configuration Ω (shown in Fig. 1c), associated with the deformed vessel representing the coupled state under fluid flow and pressure, is defined as where the deformation determines the (unknown) inner radius (r in ), the outer radius (r out ), and the vessel length (l).Finally, isochoric deformation arising from a state combined inflation, axial extension, and torsion within an elastic tube is denoted by where k = 2π 2π−α , λ z is the (constant) axial stretch and Φ is the twist angle.Note that, in ex vivo studies with coupled flow and deformation, the parameter λ z should be set based on the observed ratio of the axial length of a vessel segment before and after excision, noting that λ z > 1 due to residual stresses in vivo.For the data used in this study, λ z = 1.4 for both the control and hypertensive animals, as reported in [32].The ex vivo measurements of deformation, after introducing fluid flow through the vessel, are performed in vessels stretched and mounted to match this measured ratio.

Two-layer hyperelastic model
Within the HGO framework, a two-layer hyperelastic wall model accounting for the media (γ = M ) and adventitia (γ = A) (Fig. 1a) is formulated by representing the Cauchy stress σ = σ M + σ A as the sum of the stress in each layer [8], where Ψ γ is the Helmholtz free energy for each layer, and has the form, In (5), c γ represent the elastic moduli for the isotropic constituents (mostly elastin) in each layer, J = det(F), where F is the deformation gradient of (4), b = FF T , and Ĩlγ = A lγ : C where C = J −2/3 C (C = F T F) and A lγ = a 0lγ ⊗ a 0lγ (l = 1, 2, γ = A, M ).In (6), k 1γ and k 2γ are elastic parameters for the anisotropic constituents (mostly collagen) in each layer (Fig. 1a).Lastly, Eulerian and Lagrangian vectors, a lγ and a 0lγ (respectively), associated with collagen fiber directions are determined via (γ = A, M ), where β γ are the fixed collagen fiber angles in each layer (Fig. 1a).

Pressure-area relation
We obtain a hyperelastic pressure-area relation by integrating the radial component of the stress equilibrium equation.Neglecting inertial terms and assuming a quasi-static state this stress equilibrium equation, expressed in the current configuration, is given by where r in = r(R in ) and r out = r(R in +H), and σ rr , σ θθ are the radial and circumferential normal stress components, respectively.Here, H denotes the undeformed vessel wall thickness (Fig. 1b).Balance of forces between the transmural blood pressure and the radial component of the normal stress in the wall is enforced by the condition, Equations (4-5) are used to formulate the integrand in (9), which is evaluated with the aid of symbolic computation software (MAPLE 2019).
The resulting pressure-area relation can be written as and H M is the (reference) thickness of the media.Recall that the relation r in = a(z, t)/π is used to express the inner radius (10) in terms of the vessel area, ultimately for incorporation into (1).For brevity, the mathematical forms of the integrands F M and F A are not included here as these are lengthy expressions imported from MAPLE into MATLAB (R2021b).The integral is evaluated numerically using the MATLAB "integral" command which employs global adaptive quadrature [30] (see note in Appendix, §5.2).
This final pressure-area relation (10) contains 16 model parameters listed with units and values in the Appendix (Table 4).For a fixed set of these parameters, the model prediction of wall thickness is evaluated using equations ( 10) and ( 4) via the difference r(R in + H) − r(R in ).

Ex vivo murine data
The model developed above is calibrated to murine data made available by Chesler.
The majority of data along with detailed descriptions of the experiments can be found in the study by Tabima and Chesler [32].Data measuring lumen area and wall thickness changes with increasing transmural blood pressure were measured under ex vivo biomechanical testing in excised left pulmonary artery (LPA) vessel segments from male C57BL6 mice under control (CTL) and 10-day hypoxia-induced hypertensive (HPH) conditions [32].In both the control (CTL) and hypertensive (HPH) vessel segments, 11 measurements (i = 1, . . ., 11) are included for the relationship between vessel outer diameter (D data i ) and pressure (p data i ), and 3 measurements (j = 1, 2, 3) for the relationship between vessel wall thickness (T data j ) and pressure (p data j ).For each group, these measurements represent average values over 4 control (CTL) and 5 hypertensive (HPH) animals under controlled flow conditions with pressures in the range of 0-50 mmHg.Specific pressure values for each group are noted in Fig. 2a-b.

Fixed model parameters
Several of our model parameters can be fixed at representative values based on literature values or details of the experiments used to calibrate the models.First, we assume that the vessels have no twist, i.e.Φ = 0 • and that the opening angle in the stress-free reference state is α = 94.2 • .This latter value was obtained from literature reporting measurements in rings extracted from murine LPA vessels [37].Moreover, to mimic the in vivo setting, excised vessels were stretched to 140% of their length after extraction prior to mechanical testing [32], corresponding to a fixed value of λ z = 1.4 in (4).Since detailed histology for the murine LPA is not available, we use literature values reported in the HGO model [8] to set collagen fiber angles at β M = 29 • and β A = 62 • .Finally, we assume that the media occupies 2/3 of the vessel wall thickness in the stress-free reference state.These fixed parameter values are summarized in the Appendix (Table 4).Accounting for these six assumptions, for the parameter dependency R out = R in + H and observing that the model is independent of L yields the following 8 parameters to be estimated 2.6 Parameter estimation, sensitivity, identifiability and model reduction In the context of our model and data, we formulate a parameter estimation problem determining m parameters q * that minimize the least squares cost J as , where: where the 14-component residual vector s(q) is given by with i = 1, . . ., n 1 and j = 1, . . ., n 2 , with n 1 = 11 and n 2 = 3.The optimization problem is solved using a Nelder-Mead direct search simplex algorithm [17] minimizing J (q) in ( 13) using the routine "fminsearch" in Matlab.
Note, the mathematical model is used to evaluate the term a data in,i by first converting the outer diameter data (D data i ) to an inner radius using (4) and then using (10) to determine the values p(a data in,i ).The term T (p data j ) is evaluated as outlined at the end of §2.3.Calibration of the model to data is done in an iterative manner, gradually reducing the model complexity and number of parameters estimated using sensitivity analysis and subset selection.
Sensitivity analysis is performed after parameter estimation (with n data points) using local methods calculating the n × m sensitivity matrix χ = ∇ q s(q) using a first order finite difference scheme.Prior to calculation of sensitivity derivatives, a linear mapping is used to normalize across scales due to the diverse set of parameters and units in our model.Specifically, a perturbed interval about the kth parameter estimate ).This yields, via the Chain rule, the derivative trans- ∂y dy dx , resulting in a multiplying factor of 2αq * k in transforming raw sensitivities to their scaled counterparts.The value α = 0.1 was prescribed and all sensitivity derivatives above were approximated using first-order finite difference approximations with a step size chosen sufficiently small (10 −7 ).This choice ensured numerical convergence of all scaled parameter sensitivity derivative computations across all cases considered in this study.
Subset selection and model reduction is performed using the eigenvalue method [3,4,15,23] that analyzes the magnitude of eigenvalues and corresponding eigenvectors for the m × m Fisher information matrix approximated as χ T χ at q = q * [28].In our study, the eigenvalue subset selection method is guided by physical properties of our model at each stage of the overall process.At q = q * , our subset selection analysis and model reduction is summarized by the following iterative procedure.
1. Determine the eigenvalues of the matrix χ T χ. 2. Check if the smallest eigenvalue of χ T χ is below a specified threshold η (see, e.g.Fig. 4).3.If step 2 is satisfied, examine the eigenvector corresponding to the smallest eigenvalue.4. Mark the order 1 components of the eigenvector in step 3. 5. Parameters corresponding to the marked vector components in step 4 are potentially unidentifiable and considered as candidates for fixing at nominal values, or uncovering parameter dependencies.During the course of this iterative procedure, we ensure that the cost J (q * ) is preserved.This approach strikes a balance between model reduction and robust optimization, preserving the quality of curve-fits within the context of the given data set as the process advances.

Results
We apply the iterative approach for estimating the non-fixed parameters in (12) as outlined below.§3.1 estimates the 8 non-fixed parameters for the control animals.Results of parameter estimation, sensitivity analysis, and subset selection yield a reduced model with 6 parameters.§3.2 estimates these 6 parameters for the control animals using the reduced model.
Results of this analysis enabled further model reduction, yielding a model with equal elastic moduli in the two layers.The resulting model has 5 parameters.§3.3 estimates these 5 parameters in the further reduced model for both control and hypertensive animals.Results reveal that two parameters are correlated.Fixing one of the correlated parameters yields the final 4-parameter model.§3.4 examines the final, fully reduced 4-parameter model.To investigate effects of fixing one of the two correlated parameters identified in §3.3, we conducted a study evaluating impacts of varying the fixed parameter.We report parameter ranges for successful results, preserving quality of curve-fits via bounds on the least squares error for both control and hypertensive animals.

Baseline control animal model (8 parameters)
We first estimate the 8 non-fixed parameters for the control animals Fig. 3 Normalized parameter sensitivities for the control (CTL) animals with 8 estimated parameters across the 14 data points.
Initial and estimated parameter values for this case are reported in Table 1 (initial values are also shown Table 4 in the Appendix).
The initial value of the reference wall thickness H is set based on the first data measurement in the experiments.The initial value of R in (1 mm) is set using an order of magnitude estimate for the LPA.Setting initial values for the remaining 6 parameters is challenging given that experiments yielding direct measurements for their (nominal) values do not exist.Consequently, isotropic elastic moduli c M , c A are set to initial values of 10 kPa, an accurate order of magnitude estimate for this type of biological soft tissue.Due to nonlinearity in the anisotropic portions of the model, initial values for the remaining parameters are determined by systematic variation of initial parameter choices.Initial values for these parameters yielding a high cost J were rejected to arrive at at a combination of initial values with a curve fit to the data of good quality.The resulting combination of initial values for k 1M and k 1A is 1 kPa and 0.3 kPa, respectively, while the initial values for k 2M and k 2A are based on those reported in the HGO study [8].This combination of initial values yields the most consistent set of results across all cases considered.Model predictions with estimated parameters depicting pressure vs. area and the wall thickness vs. pressure (shown in Fig. 2a and b) provide excellent fits to the control animals.Inspection of estimated parameters reveal that the adventitia parameter value k 2A ≈ 0 (see Table 1), suggesting that we can eliminate the anisotropic terms for the adventitia (the last two terms in (5)).This observation suggests that the parameter k 1A is structurally unidentifiable since this parameter only appears in the last two (anisotropic) terms of (6) in the adventitia (γ = A).
The Fisher information matrix χ T χ is used to evaluate the eigenvalues depicted in Fig. 4a.Examination of the eigenvector of χ T χ (Fig. 4c) corresponding to its smallest eigenvalue with η < 10 −10 flags the parameter k 2A (has an order 1 component), indicating that this parameter is unidentifiable.This designation is consistent with result of sensitivity analysis (shown in Fig. 3), which demonstrates that the sensitivities for both for k 1A and k 2A are small relative to the other parameters.
Taken together, these findings motivate a model reduction in which k 1A ≈ k 2A ≈ 0. Thus, in the next step we analyze a 6-parameter reduced model eliminating the anisotropic terms for the adventitia in the stress-strain law.

Reduced control animal model (6 parameters)
Parameter values are initialized as described in §3.1.The 6 parameters to be estimated for the control animals are Results are listed in Table 1.Again, for the control animals the quality of curve fits of the model to the pressure vs. area data (Fig. 5a) and the wall thickness vs. pressure data (Fig. 5b) is preserved, with a slight reduction in overall cost from J = 1.707 • 10 −4 to J = 1.674 • 10 −4 .In addition, the estimated values are preserved within 0.25% for the geometric parameters (R in , H), concurrent with a 3.0% reduction in the elastic modulus k 1M , a 11.5% reduction in the elastic modulus c M and a substantial increase (2.4×) in the modulus c A .Finally, the dimensionless parameter k 2M increases by 1.5%.Examination of the eigenvector of χ T χ (Fig. 6h) corresponding to its smallest eigenvalue with (η < 10 −7 ) shown in Fig. 6g flags the two parameters c A and c M with order 1 components; c A is the dominant component.The sensitivity for c A is also small relative to the other parameters (Fig. 6a-f).These findings, combined with the observation that c M and c A exhibit the largest changes in estimated values between the 8-parameter and 6-parameter fits, motivate a further reduced 5-parameter model examined in the next stage of the process.
Since two of the remaining 5 parameters are geometric parameters, we retain three elastic parameters describing the isotropic and anisotropic responses in the model.Thus, the only elastic parameters estimated in the next step are c M , k 1M and k 2M .Specifically, the reduced model analyzed in the next section has elastic moduli in the media and adventitia that are set equal, whereas the two layers still retain distinct collagen fiber angles (β M , β A ).

Reduced control and hypertensive animal model (5 parameters)
In the reduced 5-parameter model, elastic parameters in the two layers are assumed equal, i.e., c A = c M , k 1A = k 1M , and k 2A = k 2M .The parameter vector estimated for this model is This model is fitted to data from both the control (CTL) and hypertensive (HPH) animals.Values of the five model parameters are initialized as described in §3.1 and the estimated parameter values are reported in For the control (CTL) animals, the quality of curve fits of the model to the pressure vs. area data (Fig. 7a) and the wall thickness vs. pressure data (Fig. 7b) are preserved, with a slight reduction in overall cost from J = 1.674 • 10 −4 to J = 1.641 • 10 −4 .The curve fit for hypertensive (HPH) animals has a significantly lower cost (J = 0.778 • 10 −4 ) due, in part, to the smaller range of variation in the pressure-area curve caused by vessel wall stiffening (Fig. 7a & Table 2).In the hypertensive animals, geometric parameters exhibit a small increase in vessel wall inner radius R in (≈ 394µm vs. ≈ 377µm) and a large (and expected) increase in reference wall thickness H (≈ 62µm vs. ≈ 45µm).The altered dynamics in the hypertensive animals are reflected by a substantial increase in the elastic modulus k 1M (13.38kPa vs. 0.60kPa) and in the dimensionless parameter k 2M (1.30 vs. 0.58), both associated with collagen stiffness.Concurrently, in the hypertensive animals there is a substantial drop in the elastic modulus c M (21.50kPa vs. 0.58kPa), associated with elastin deformation, relative to the control animals.For this model, subset selection evaluating the eigenvalues and eigenvectors of χ T χ (Fig. 7d-g) for the estimated parameters (Table 2, Fig. 7d-g) reveals that the eigenvectors for the smallest eigenvalue (Fig. 7d,f) flag two parameters k 1M and c M as the dominant components in the control (CTL) and hypertensive (HPH) animals, respectively (Fig. 7e vs. Fig.7g).These observations suggest a parameter dependency between the two elastic moduli c M and k 1M , which share the same units, that is evaluated in the next section using our final 4-parameter model.

Parameter dependencies and range estimates for the final model (4 parameters)
To study the implications of fixing one of the two correlated parameters, for both the control and hypertensive animals, our final model fixes c M while still estimating k 1M , i.e., analysis in this section estimates the 4 parameters, Since we do not have data from independent experiments for k 1M , we repeat optimization while varying this parameter.To preserve quality of the curve fits, parameter ranges were determined by enforcing the condition that the cost increased by no more than 10% (to two decimal places) relative to the values in Table 2, and across several runs of the optimization.Based on this criteria, the maximum costs used as cutoffs are set to J = 1.81 • 10 −4 (CTL) and J = 0.86 • 10 −4 (HPH).The resulting estimated parameter ranges are reported in Table 3.The corresponding curve fits are not shown as they were visually identical to those shown in Fig. 7.For both the control and hypertensive animals, the estimated parameter values for k 1M varied inversely with c M .For the hypertensive animals, the geometric parameter ranges exhibit a small increase in vessel wall inner radius R in and a large increase in reference wall thickness H.The ranges of values for the elastic moduli associated with collagen (k 1M , k 2M ) are substantially higher while the range of values for the modulus associated with elastin (c M ) is substantially lower, compared to the control animals.Furthermore, the range for c M includes zero, indicating that the hypertensive model can be fit well to the data with a negligible biomechanical contribution of elastin to the stress-strain and wall thickening response to increasing pressure.This study develops a new mechanistic wall model relating transmural blood pressure and vessel lumen area using a two-layer nonlinear hyperelastic HGO model incorporating residual stresses and anisotropy.The new model is calibrated using pressure and wall thickness data from [32] for control and hypertensive mice.
In the hypertensive animals, pulmonary hypertension is induced by placing the animals in a hyperbaric chamber exposing them to hypoxia for 10 days.Model calibration and systematic model reduction is achieved by combining sensitivity analysis, subset selection, and parameter estimation.The results demonstrate that this detailed structural continuum mechanics model, containing a large number of parameters, can be systematically reduced to capture differences in key model parameters between control and hypertensive animals.
To our knowledge, this is the first study to carry out robust parameter estimation and model reduction by simultaneously predicting the increase in lumen area and the decrease in wall thickness as pressure is increased.The results reveal that these biomechanical responses can be accurately captured using a model that retains a single set of three elastic moduli delineating the contributions of collagen and elastin under the loading protocol of the associated experiments.Specifically, the material parameter associated with elastin (c M ) was the dominant contributor to nonlinear stiffening in the control animals.By contrast, in the hypertensive animals, the contribution of c M was negligible and nonlinear stiffening is dominated by material parameters associated with collagen (k 1M , k 2M ).Taken together, these findings are consistent with well-known increases in collagen content in the wall of large pulmonary arteries with hypoxia-induced PH [20,35].Note that our analysis uses the same opening angle value (α = 94.2 • [37]) for both the control and hypertensive models, based on the only known measurements of this quantity in a similar vessel and species after 10-days of hypoxic conditions [9] (see Fig. 7 therein).
The robustness of our model and approach is evidenced by its accurate and simultaneous prediction of both pressure and wall thickness changes under deformation, for both control and hypertensive data sets.Our systematic approach to parameter identifiability, subset selection and model reduction decreased the overall number of parameters in the model while preserving the quality of curve-fits to the data at each stage of the iterative procedure.
Limitations include common parameter estimation challenges when the number of model parameters and/or variables increases relative to the variables or quantities of interest for which data is available, as well as the lack of known nominal values for some model parameters in large pulmonary arteries.One challenge is non-uniqueness of parameter estimates due to the infeasability of guaranteeing a solution of the optimization problem that is a global minimum of the cost function across the parameter landscape.A second challenge is the local nature of sensitivity measures underlying the identifiability techniques used in this study, i.e. the final reduced model is not guaranteed to be unique.Indeed, the accuracy and robustness of the approaches presented in this study can be enhanced through both extended ex vivo and in vivo studies informed by the model presented here.Our approach for sensitivity and identifiability analysis and model reduction is rooted in prior works describing parameter subset selection techniques [3,4,22,23] using an eigendecomposition of the matrix χ T χ, but similar results could likely be obtained using other methods.While global sensitivity analysis techniques exist [10,29,13], most subset selection techniques are local.The method for identifiability analysis used here is based on eigenvalues but, as discussed in several previous studies, similar results can be obtained using other methods [7,15,19].Overall, sensitivities or unidentifiable parameters for particular variables of quantities of interest can suggest which types of data will be most influential in an expanded data set.Where practical, examples of extensions include augmentation of ex vivo biomechanical testing to include measurement of the vessel opening angle, as well as incorporation of in vivo data measuring BP, flow and lumen area prior to sacrifice of the animal(s).
Our model is also based on assumptions of hyperelastic deformation and geometric idealization of the stress-free reference state (Fig. 1b) as a segment of a cylindrical ring; in reality, the vessel wall may exhibit viscoelastic effects under pressurization and/or deviate from circular arcs in the cut open rings.The new pressure-area relation developed in this study has the potential for incorporation into 1D cardiovascular network models of coupled fluid-solid dynamics in large pulmonary arteries.Some possible approaches include direct incorporation and coupling of the pressure-area relation within the 1D fluids network solver or, alternatively, using the pressure-area relation as a high fidelity model for emulation using simpler empirical models [11,18,33] or statistical models.Overall, the techniques and findings presented in this study demonstrate the potential for development and systematic reduction of more realistic models of key relations (e.g.pressure-area) through the integration of data-driven mathematical approaches for ex vivo experiments with modeling approaches predicting in vivo dynamics in cardiovascular biomechanics.

Appendix
Table 4 List of all parameters for the model developed in this study.Estimated parameters listed with an asterisk (*) are ultimately fixed (during the course of the model reduction).Note that the model no longer depends on the parameter L when the twist angle Φ is assumed to be zero (see ( 4)).

Complete set of model parameters
For convenience, the full set of model parameters, their descriptions, units, designation of parameter type (estimated, fixed, dependent or eliminated) and the associated fixed or initial values are summarized in Table 4.

Pressure-area relation
A Github link to Matlab functions for the integrands F M (r in , r) and F A (r in , r) in ( 10) is available upon request.

Fig. 1
Fig.1Foundations of the nonlinear hyperelastic wall model: (a) Illustration of a cross-section of a large artery wall (redrawn from[8]); (b) the stress-free reference state Ω 0 defined in equation (2) where R in is the inner radius, Rout is the outer radius, H is the wall thickness, L is the axial length and α is the opening angle; (c) the current configuration Ω defined in(3).Note that the (deformed) inner radius (r in ), outer radius (rout) and axial length (l) are all determined after the model equations are solved.

Fig. 2
Fig. 2 Results from estimating 8 parameters (listed in Table 1) for the control (CTL) animals.(a) pressure vs. area, model predictions of the outer area (black) vs. data (circles) and the inner area (red); (b) wall thickness vs. pressure model predictions compared to the 3 data points (squares); (c) the residual vector (14) across the 14 data points.

Fig. 4
Fig. 4 Identifiability results using the eigendecomposition of the information matrix (χ T χ) for the control (CTL) animals with 8 estimated parameters: (a) log-plot of the eigenvalues of χ T χ; (b) components of the eigenvector of χ T χ corresponding to the smallest eigenvalue of χ T χ that is less than η = 10 −10 (asterisk).

Fig. 5
Fig. 5 Parameter estimation results for the reduced model for the control (CTL) animals with 6 estimated parameters: (a) pressure vs. area model predictions of the outer area (black) vs. data (circles) and the inner area (red); (b) wall thickness vs. pressure model predictions compared to the 3 data points (squares); (c) plot of the residual vector (14) across the 14 data points.

Fig. 6
Fig.6Identifiability results computed using eigendecomposition of the information matrix (χ T χ) for the reduced 6 parameter model with data from the control (CTR) animals: (a-f) normalized parameter sensitivities for the 6 estimated parameters across the 14 data points; (g) log-plot of the eigenvalues of χ T χ; (h) components of the eigenvector of χ T χ corresponding to the smallest eigenvalue of χ T χ less than η = 10 −7 (asterisk).

Fig. 7
Fig. 7 Parameter estimation and identifiability results for the reduced 5-parameter model for the control (CTL) and hypertensive (HPH) animals: (a) pressure vs. area model predictions of the outer area (black) vs. data (circles) and the inner area (red); (b) wall thickness vs. pressure model predictions compared to the 3 data points (squares); (c) plot of the residual vector (14) across the 14 data points; (d,f) log-plot of the eigenvalues of χ T χ in the normotensive (d) and hypertensive (f) animals; (e,g) components of the eigenvector of χ T χ corresponding to the smallest eigenvalue of χ T χ (asterisk) in the control (e) and hypertensive (g) animals.

Table 1
Estimated parameter values for the control (CTL) animals with the 8-parameter model (column 5) and the reduced 6-parameter model (column 6).

Table 2 .
For comparison, results of the 6-parameter model are also included in the table.Note that the estimated parameter values of c M , k 1M , and k 2M for this reduced model should be interpreted as aggregate elastic parameters for the entire vessel wall, i.e., representing all layers.

Table 2
Estimated parameter values for the reduced 5-parameter model for control (CTL) (column 6) and hypertensive (HPH) (column 5) animals.For comparison, results for the reduced 6-parameter model are also shown (column 4).

Table 3
Estimated parameter ranges for the final model in both the control (CTL) and hypertensive (HPH) animals based on optimization with 4 parameters.The cost J was allowed to increase by no more than 10% in establishing the estimated parameter ranges.