Uncertainty quantification in subject-specific estimation of local vessel mechanical properties

Quantitative estimation of local mechanical properties remains critically important in the ongoing effort to elucidate how blood vessels establish, maintain, or lose mechanical homeostasis. Recent advances based on panoramic digital image correlation (pDIC) have made high-fidelity 3D reconstructions of small-animal (e.g., murine) vessels possible when imaged in a variety of quasi-statically loaded configurations. While we have previously developed and validated inverse modeling approaches to translate pDIC-measured surface deformations into biomechanical metrics of interest, our workflow did not heretofore include a methodology to quantify uncertainties associated with local point estimates of mechanical properties. This limitation has compromised our ability to infer biomechanical properties on a subject-specific basis, such as whether stiffness differs significantly between multiple material locations on the same vessel or whether stiffness differs significantly between multiple vessels at a corresponding material location. In the present study, we have integrated a novel uncertainty quantification and propagation pipeline within our inverse modeling approach, relying on empirical and analytic Bayesian techniques. To demonstrate the approach, we present illustrative results for the ascending thoracic aorta from three mouse models, quantifying uncertainties in constitutive model parameters as well as circumferential and axial tangent stiffness. Our extended workflow not only allows parameter uncertainties to be systematically reported, but also facilitates both subject-specific and group-level statistical analyses of the mechanics of the vessel wall.


INTRODUCTION
Detailed quantification of the nonlinear, anisotropic mechanical behavior exhibited by blood vessels is essential for understanding the wall mechanics, for informing models of the associated mechanobiology, and for building appropriate fluid-solidinteraction models of the hemodynamics. Such detail is generally possible only for excised specimens subjected to well-defined boundary conditions. Whereas quantification is ultimately desired for human arteries, mouse models of vascular development, health, aging, and disease continue to prove increasingly valuable, particularly given the ability to collect longitudinal data on wild-type, genetically modified, pharmacologically treated, and surgically manipulated mice in large numbers. Given that many normal segments of elastic and muscular arteries are nearly cylindrical, miniaturized biaxial pressure-distension and axial force-extension testing devices have provided considerable data for such quantification. [1][2][3] There are many cases, however, wherein the geometry is complex, such as in aneurysms, atherosclerosis, and dissections, which necessitates advanced experimental methods. We thus developed a panoramic digital image correlation (pDIC) method than can be coupled with optical coherence tomographic (OCT) examinations of murine vessels to quantify material properties regionally. 4,5 This approach has since provided increased insight into diverse conditions, including thoracic aortic aneurysm, abdominal aortic dissection, and aortic tortuosity. [6][7][8] Our pDIC approach yields high-fidelity 3D reconstructions of murine vessels imaged in multiple quasi-statically loaded configurations, which allows computation of the full field of local deformations along the vessel surface. 4,9 To estimate associated best-fit values of constitutive model parameters regionally, we use a novel inverse characterization approach based on the principle of virtual power. 10 Adding local thickness measurements derived from OCT allows the inverse characterization to be extended from DIC to digital volume correlation, as illustrated in our study of aortic dissections that include significant intramural thrombosis. 6 Notwithstanding these prior technical advances, there has remained a pressing need for more robust quantification of biomechanical properties, complete with estimates of their associated uncertainties, to enable reliable comparisons of results both across complex lesions and across mouse models. The lack of quantitative uncertainty estimates have made it unfeasible to address questions regarding inferred biomechanical properties on a subject-specific basis, such as whether stiffness differs significantly between multiple material locations on the same vessel, or whether stiffness differs significantly between multiple vessels at a corresponding material location. In previous studies, [5][6][7][8] we have circumvented this issue by performing only group-level analyses of pooled regional results, while also limiting our analysis to material points at which the constitutive model explained a very high proportion of the variance in the observed mechanical data (as quantified by the coefficient of determination). We note, however, that the ability of the model to fit the data well does not necessarily correspond to a high degree of identifiability in the material parameter values, or similarly in derived estimates of the vessel wall's mechanical properties, such as tangent stiffness or stored energy. 3 Hence, there was a critical need to extend the modeling framework to provide measures of uncertainty that can be reported alongside the locally estimated material parameters and related quantities of interest.
In the present study, we have integrated a novel uncertainty quantification and propagation pipeline within our inverse modeling approach, relying on empirical and analytic Bayesian techniques. To demonstrate the approach, we present illustrative results for three murine ascending thoracic aorta specimens, quantifying uncertainties in constitutive model parameters as well as circumferential and axial tangent stiffness components (Section 3). We note that consistent comparisons of data across different mouse models promises to provide far greater insight into vascular health and disease (cf. Bellini et al. (2017), 11 Humphrey & Tellides (2019) 12 ) than detailed studies that only compare wild-type to a single model of interest, as is common. Toward this end, herein we applied our pipeline to data from three different mouse models: wild-type control, apoliproprotein E-null (Apoe −∕− ) infused with angiotensin II (Ang II), and wild-type treated with -aminopropionitrile (BAPN). In addition to facilitating both subject-specific and group-level statistical analyses of the mechanics of the vessel wall, our extended workflow will allow mechanobiological models that utilize estimated tissue-level mechanical properties downstream of the present pipeline to account for their associated uncertainties.

Summary of approach
Herein, we augment previously established procedures for mechanical testing, high-resolution imaging, and inverse modeling with a new integrated workflow for uncertainty quantification and propagation ( Figure 1). As before, we first collect standard biaxial mechanical data to examine the gross response of the vessel following established protocols, 3 which enables estimation of the complete reference configuration (Section 2.2.2). A speckle pattern is then applied to the vessel, which is mounted within a custom device for panoramic imaging. The acquired images are post-processed using pDIC, which yields 3D reconstructions of the vessel in each imaged configuration as well as the complete field of local deformations. Combining the deformation and loading data with local estimates of the wall thickness derived from OCT (Section 2.3) allows inverse modeling to estimate local material model parameters (Section 2.4), from which we can compute mechanobiologically relevant quantities of interest, such as tangent stiffness. To advance this approach, we now quantify uncertainties in the locally estimated model parameters using a Bayesian framework (Section 2.4.1), and then propagate those uncertainties to obtain posterior distributions of our quantities of interest at each point along the vessel wall (Section 2.4.2).

FIGURE 1
Summary of the complete workflow required to estimate mechanical properties locally. Boxes shown in blue highlight contributions of the present study, which focuses on the quantification of uncertainty in material model parameters as well as the propagation of that uncertainty to broader quantities of interest (e.g., circumferential and axial stiffness) whose estimated values depend on multiple model parameters. config., configuration; OCT, optical coherence tomography; pDIC, panoramic digital image correlation

Illustrative data set
To demonstrate and validate our novel inverse modeling and uncertainty quantification pipeline, we applied our approach to three ascending thoracic aorta specimens excised from male mice with a C57/BL6 genetic background. To highlight potential differences in the uncertainty of estimated mechanical properties across different specimen types, we herein show results for one wild-type control specimen, one Apoe −∕− specimen infused with Ang II for 14 days to induce an aneurysm coincident with hypertension, 13 and one juvenile wild-type specimen treated with BAPN for 14 days to induce a dissection caused by inhibition of collagen cross-linking (Table 1). 14 All animal protocols were approved by the Yale University Institutional Animal Care and Use Committee and followed methods detailed previously. 3,15

Biaxial mechanical testing and panoramic imaging
Noting the crucial importance of biaxial data for the accurate estimation of vessel mechanical properties, we collected gross pressure-diameter and axial force-length data for each vessel following established procedures. 3 Briefly, each vessel was cannulated, secured with sutures, and mounted within a previously validated computer-controlled testing system that permitted both cyclic distension and axial extension. 16 Following preconditioning, the vessel diameter was measured at pressures ranging from 10 to 140 mmHg while holding the axial stretch constant at 0.95, 1, and 1.05 times the initially estimated in vivo stretch. Similarly, the axial length was measured at axial forces ranging from 0 to 35 mN while holding pressure constant at 10, 60, 100, and 140 mmHg. In preparation for subsequent mechanical data collection within the panoramic imaging system, the estimated stretch value at which axial force would vary the least during pressurization was adopted as an updated estimate of the in vivo stretch ( iv ). 17 Each vessel was air-brushed with black and white India ink to generate a random speckle pattern over the entire outer surface, then mounted within a panoramic imaging system previously described. 4,9 Panoramic images were acquired while each vessel was loaded quasi-statically to known values of axial stretch (0.95 iv , iv , 1.05 iv ) and pressure (10-140 mmHg, in steps of 10 mmHg), totaling 42 unique configurations ( Figure 2). In each configuration, vessels were imaged from eight views, thus yielding 336 images total per vessel. Due to the high number of images acquired, post-processing of the imaging data for subsequent modeling purposes has historically required about eight hours of interactive input per specimen. In anticipation of the present study and others, we systematically automated this post-processing pipeline, reducing the total time required to less than 30 minutes per specimen with no change in the quality of the results. These automation efforts are described in the Appendix.
To arrive at 3D reconstructions of each vessel in every configuration, we performed pDIC between pairs of post-processed images, following previously described methods. 4,9 Note that the set of 3D reconstructions for each vessel have material point correspondence, enabling the computation of the complete local deformation (tensor) field at every point on the adventitial surface, which is crucial for subsequent inverse modeling and uncertainty quantification ( Figure 2; Section 2.4). Each vessel surface was discretized using a 40 × 25 array of quadrilateral elements as in Weiss et al. (2020) 8 prior to inverse modeling.

Local vessel wall thickness estimation
Accurate (local) material model parameter estimation is highly dependent on knowledge of (local) specimen thickness, which is inversely proportional to the estimated stress values used during regression. We recently integrated OCT-derived local thickness measurements within our pDIC-based inverse modeling pipeline, allowing us to no longer rely on simplifying uniform-thickness assumptions. 6 In the present study, we adopt the approach described in Bersi et al. (2019), 6 but make several practical improvements to the processing of the data as described below. OCT data were collected for each vessel in the reference configuration  6 to form a complete cross section orthogonal to the central axis ( Figure 3).
As an improvement upon previous approaches that required interactive tracing of the inner and outer vessel wall at each cross section, we developed an automated algorithm to extract these contours, thus dramatically reducing the post-processing time as well as eliminating variability resulting from observer subjectivity. First, the image is thresholded using Otsu's method, 18 such that high-intensity pixels are categorized as likely tissue regions and low-intensity pixels as likely non-tissue regions. Second, the image data is interpolated on a polar coordinate system, whose origin is at the mean location of all high-intensity pixels ( Figure 3).
The medial contour of the vessel wall is then obtained via a locally estimated scatterplot smoothing (LOESS) regression of the high-intensity pixel coordinates in polar space. 19,20 Periodicity was enforced by performing the regression on data which were repeated in the circumferential direction. The local regression model was chosen to be a second-order polynomial with a tricube weight function and span of 90 • in the circumferential direction. To remove the outlier effects of spurious bright pixels far from the true vessel wall, only points with radius between 1/3 and 3 times the median radius were included in the regression. The regression was further rendered robust to outliers by performing a second iteration in which the data were weighted using a bisquare function in the radial direction, with a span of 6 times the median absolute radial deviation in the first iteration. 19 At each circumferential coordinate value , we computed the minimum and maximum radius residuals, min, and max, , from the medial contour. To extract the inner and outer contours, we performed the same LOESS procedure as before, but on the min, and max, respectively, and added the result to the medial contour. The local thickness was finally computed as the difference

FIGURE 3
Estimation of local thickness via optical coherence tomography (OCT). Images are acquired at 100 cross-sections (only 20 highlighted above for clarity) along the length of the vessel in its (loaded) reference configuration, and each image is interpolated on a polar coordinate system with the origin at the mean location of all pixels with intensity above a minimum threshold. The medial (magenta), inner (blue), and outer (red) contours of the vessel wall are then estimated using local polynomial regression (LOESS) of high-intensity pixels with radius greater than a minimum value (green) between the inner and outer radius, using the approximation that the surface normal nearly coincides with the radial coordinate direction at all points (this approximation holds except, for example, in the presence of large saccular aneurysms). After the contours of all cross sections are obtained, the 3D outer contour stack is registered to the pDIC-derived reconstruction of the reference configuration using a least-squares-optimal rigid body transformation as described in

Bayesian parameter estimation
In line with several of our previous studies, 3,5,6,21 we modeled the vessel wall as a four-fiber-family hyperelastic constrained mixture of (1) isotropic elastin, (2) smooth muscle cells (SMCs) and circumferential collagen, (3) axial collagen, and (4) diagonal collagen, with strain energy density of the form where is the right Cauchy-Green deformation tensor; It is generally assumed that all of the constituent mass fractions can be accurately prescribed from histology. 5,6 Therefore, the vector of model parameters remaining to be estimated at each material point is which are all strictly positive. Adopting a Bayesian framework, we may state our model parameter estimation objective as follows: given the observed mechanical data, and taking into account prior knowledge regarding the model parameter values, where is the maximum a posteriori estimate of , ( | exp ) is the posterior PDF of given exp , ( exp | ) is the likelihood of observing the data given , and ( ) is the prior PDF of . (Note the shorthand notation for brevity to denote any PDF, though formally these would be symbolically distinguished since they do not have the same form in general. In this shorthand notation, the form of each should thus be inferred from the function argument. In some cases, for clarity, we will distinguish PDFs of different forms using subscripts on .) Equivalently, We assume for analytic convenience, as in a least-squares regression, that the residuals exp, − th, are independent and follow a zero-mean multivariate normal distribution with some covariance matrix, which we treat as a nuisance parameter.
For simplicity, the prior distribution of can also be modeled as a multivariate normal distribution. When there is no prior knowledge, a flat prior may be chosen ( ( ) ∝ 1); in this case, the maximum a posteriori estimate coincides with the maximum likelihood estimate. Equivalently, the log-prior is chosen to be either quadratic with respect to or constant, depending on the FIGURE 4 Near a maximum a posteriori value , the logarithm of the posterior probability density function (PDF)-or equivalently, the log-likelihood function alone if a flat prior is used-can be approximated well by a quadratic function ("Laplace's approximation"), corresponding to a normal approximation of the posterior distribution. [22][23][24] For one (multiple) parameter(s), the negative inverse of the curvature (Hessian matrix) corresponds directly to the variance (covariance matrix) of the approximated posterior distribution availability of prior knowledge. Under these assumptions, it can be shown that at = , the sum of the log-likelihood and logprior (equal to the log-posterior ln( ( | exp )) up to an additive constant) is also asymptotically quadratic with respect to ; this is known as "Laplace's approximation." [22][23][24] Since the logarithm of the normal PDF is quadratic, it follows that the posterior distribution of is asymptotically normal, with mean and covariance matrix = − −1 , where  the Hessian matrix of the log-posterior ( Figure 4).
Since the posterior distribution of the model log-parameters is multivariate normal, the marginal distributions of each are also normal, with mean and standard deviation = √ Σ . The marginal posterior distributions of the untransformed model parameters Ξ are thus log-normal, with median and geometric mean g, = , and geometric standard deviation (GSD) g, = , which can be used as a scalar summary of the relative uncertainty in a given model parameter. Note that > 0 ⇐⇒ g, > 1, with g, → 1 denoting complete certainty. From these geometric distribution parameters, equi-tailed credible intervals for Ξ can be constructed of the form [ g, ∕ g, , g, g, ], where defines the coverage probability (for a 95% credible interval, ≈ 1.96).
Alternatively, highest posterior density credible intervals can be constructed numerically by searching for the narrowest interval that yields the desired coverage probability. 25 Using a flat prior, we applied the above inverse modeling framework pointwise along the discretized vessel surface, yielding a set of location-specific posterior distributions for the estimated set of model log-parameters. Due to the cost of computing the gradient of the objective function, we started the optimization by performing 10 iterations using a derivative-free genetic algorithm, 26 and then used the best current candidate to initialize a local maximizer using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, 27,28 which executed until convergence. The Hessian matrix  was computed at the optimum using finite differences.

Uncertainty propagation
Due to the phenomenological nature of the constitutive model, physical interpretability of the specific parameter values is limited, except perhaps in the case of deposition stretches and the angle of the diagonal collagen fiber family. Rather, quantities derived from the evaluation of the constitutive model (conditional on a set of parameter values) are often of primary interest, as they can be related directly to mechanobiological phenomena. 29,30 These quantities include direction-specific material stiffness, stress components, and stored energy, among others. 8,31 In addition to uncertainties in the model parameters themselves, it is thus crucial to report uncertainties in these derived quantities of interest, especially for the purposes of model selection, future experimental design, computational (e.g., finite element) simulations, and application studies involving subject-specific and/or tissue-location-specific hypothesis testing. Extending the parameter estimation and uncertainty quantification pipeline where E[⋅] denotes the expected value and tr(⋅) denotes the trace. Note that T T = ∑ . By Isserlis' theorem, Therefore (exploiting the symmetry of and ), In the case of a vector function  ( ), the mean and variance of each component can be independently computed following Equations 6 and 8, while the covariance between components can be similarly derived to yield Let ln(( 1111 ) , ) denote the circumferential log-stiffness at material point for specimen (note that we again use a logarithmic transformation to map the strictly positive quantity ( 1111 ) , to an unbounded space). At the group level, we assume for simplicity that the distribution of log-stiffness is normal, with mean group and variance 2 group . For certain suitable priors on [ group , 2 group ], if the true subject-specific ln(( 1111 ) , ) were known exactly, the posterior PDF of the group-level distribution parameters would have a straightforward analytic form. For example, under a normal-inverse-gamma prior, which is conjugate to the normal distribution, the posterior distribution of [ group , 2 group ] would also be normal-inverse-gamma. 24 The same is true for the uninformative prior proportional to 1∕ 2 group , albeit with different distribution parameters. The marginal distribution of group would be a non-standardized Student's distribution (̂ ,̂ 2 , ) in both of these cases. If the uninformative prior is where is a dummy variable for integration over the entire domain of ln(( 1111 ) , ), ( group | ) is the PDF of the nonstandardized Student's distribution (̂ ( ),̂ 2 ( ), ), and  ( ) is the PDF of the -variate normal distribution  (̃ ,̃ ).

The integral in Equation 10
has no closed-form expression, but is straightforward to compute numerically by quadrature or through Monte Carlo sampling of  (̃ ,̃ ), where the result is simply a weighted average of Student's PDFs. Note that the same approach can be used to compute the bivariate posterior distribution of [ group , 2 group ], whose PDF ( group , 2 group ) is a weighted average of normal-inverse-gamma PDFs.
While we have outlined the above approach for an individual material point , we note that the analysis could be broadened to an anatomically relevant region of the vessel wall (e.g., anterior, dorsal, etc.). Similarly, for a future specimen * , the PDF for the predictive distribution of the subject-specific estimate ln(( 1111 ) , * ) can be computed by marginalizing out both group and where  (ln(( 1111 ) , * ) | group , 2 group ) is the PDF of the normal distribution  ( group , 2 group ).

FIGURE 5 Geometric standard deviation (GSD) maps of material parameters of elastin ( e ), circumferential smooth muscle cells (SMCs)/collagen ( m ), axial collagen ( a ), and diagonal collagen ( d ) for an illustrative wild-type specimen
For both the group-level mean and the new-sample predictive distribution, which are both defined above for the log-stiffness, the PDFs of the corresponding untransformed stiffness values are computed using standard distribution mapping techniques. In the general case, if = , where is a random "log"-quantity distributed according to a PDF ( ) ( ), the PDF of is which has strictly positive support.

Uncertainties in material parameters
Application of our novel inverse modeling and uncertainty quantification pipeline to ascending thoracic aorta specimens demonstrated that the resulting uncertainties are heterogeneous across different model parameters as well as throughout different regions of the vessel wall ( Figure 5). Notably, uncertainties associated with the material parameters of elastin ( e ) and diagonal collagen ( d ) were substantially higher than those associated with circumferential SMCs/collagen ( m ) and axial collagen ( a ), suggesting that the properties of elastin and diagonal collagen are more difficult to quantify using the examined range of pressure and axial stretch. Nevertheless, GSD values were typically between 1.05 and 1.4, suggesting that in most cases the unobserved/"true" parameter values are close to the maximum a posteriori estimate. In the case of a relatively high GSD of 1.4, construction of commonly reported credible intervals indicated that the unobserved material parameter value falls between ∼0.5 and 2 times its posterior median with at least 95% probability ( Figure 6). This trend indicates that overall stiffness may be estimated with high credibility even when the model parameters themselves are difficult to estimate; broadly speaking, this is because the quality of the constitutive fit to the mechanical data is influenced jointly by all of the model parameters, whose individual contributions are sometimes weak or partly redundant.

FIGURE 7
Geometric standard deviation (GSD) maps of circumferential and axial stiffness for an illustrative wild-type specimen, determined by propagating uncertainties in the model parameters through either Monte Carlo sampling or a Taylor series expansion. Under the observed uncertainties, the two methods yield comparable results. Note also that the uncertainties in overall stiffness are often substantially smaller than uncertainties in the individual material parameters (compare to Figure 5, which has different color axis limits)

Comparison between specimen types
We applied our pipeline successfully to both a wild-type control specimen as well as two non-control specimens-one from an Apoe-deficient mouse infused with Ang II and one from a wild-type juvenile mouse treated with BAPN ( Table 1). The resulting uncertainties associated with circumferential and axial stiffness components varied substantially across specimen type (Figure 8).
In general, uncertainties were smallest for the wild-type control specimen, indicating that the model parameters and derived stiffness quantities were more identifiable, despite the constitutive model's general ability to fit non-control mechanical data well. 6 The BAPN specimen, which experienced less circumferential stretch during mechanical testing due to its prior dilatation,  36 reinforced the finding that uncertainties in axial stiffness were consistently larger than those in circumferential stiffness (Figure 8b).

Uncertainties at the group level
To demonstrate the propagation of uncertainties to the level of group statistics, we computed the posterior distribution of the group-wide geometric mean of the circumferential stiffness ( g (( 1111 ) )) at a materially correspondent point across all three specimens (Figure 9a). Note that we present these results merely to illustrate our uncertainty propagation pipeline, without implying that specimens of these different types would/should be pooled in an application-focused study. Following the assumptions and methods described in Section 2.4.2, and starting from the subject-specific posterior distributions of ( 1111 ) , (Figure 9b), the resulting posterior PDF of g (( 1111 ) ) is a weighted mixture of non-standardized Student's PDFs in the log-space mapped to the linear space using Equation 12 (Figure 9c). Similarly, we computed the predictive distribution of new subject-specific estimates of ( 1111 ) , * , which is substantially more variable than the group mean (Figure 9d). For the selected point , we found that the shape and effective support of this predictive distribution was similar to the distributions of subject-specific circumferential stiffness estimates over the entire vessel wall (Figure 9e). Broadly speaking, this result indicates that point is not expected to be substantially different from the vessel as a whole.

DISCUSSION
There is now an extensive literature on nonlinear material properties in health and diverse diseases for arteries from both humans and myriad animal models. 37,38 In most studies, the associated data and quantification focus on bulk, or homogenized, properties, not site-specific variations. In many cases, however, the nidus of disease is highly localized, with focal differences in material properties nucleating atherosclerotic lesion development, aneurysmal expansion, dissection, and so forth. There is, therefore, a pressing need for more localized descriptions of arterial properties, both within the same lesion and across different examples thereof. Motivated thus, we proposed a pDIC + OCT-based approach for quantifying material properties locally in excised samples subject to in vivo relevant mechanical loading. To exploit such information, however, there remained a pressing need both to streamline and automate the data analysis pipeline and to improve the methods of quantification, particularly to assess uncertainty.
In the present study, we have augmented our pDIC + OCT approach to quantify regional material properties with a novel uncertainty quantification and propagation pipeline. We illustrated the utility of this pipeline for ascending thoracic aorta specimens from both a wild-type control mouse as well as two disease models to demonstrate its practical use, and showed illustrative results both for the uncertainties of individual material parameters and the propagation of those uncertainties to tangent stiffness at the subject-specific and group levels. This development represents an important extension of our experimental-computational approach and follows previous milestones in the establishment of our current workflow, including the use of eight panoramic views to ensure high accuracy in pDIC-based reconstructions for strain quantification, 4 the development of an inverse modeling approach based on the virtual fields method to infer point-wise material parameters, 5 and the integration of OCT-derived local thickness measurements. 6 The contribution of the present study makes it possible, for the first time, to assess the identifiability of material parameters (Section 3.1) and the uncertainties in any derived quantities of interest, such as tangent stiffness (Section 3.2), locally and on a subject-specific basis. The opportunity to account for subject-specific uncertainties also improves the credibility of group-level statistical analyses (Section 3.4).
While the presented experimental and analytic methods provide a valuable estimate of subject-specific mechanical properties, and especially their spatial heterogeneity as well as associated uncertainties, several practical limitations remain. Most notably, the mechanical testing and image acquisition techniques rely on excised samples, thus precluding longitudinal estimation of mechanical properties for the same subject. While methods to reconstruct the geometry and estimate local deformations of cardiovascular tissues from in vivo images are emerging, 7,39,40 the spatial and temporal resolution provided by current technology is not yet sufficient to eliminate the need for ex vivo pDIC-and OCT-based measurements in murine samples. Moreover, estimation of local material parameters directly from in vivo images remains difficult, even in human subjects which have a considerably larger vessel geometries. [41][42][43][44][45] An additional caveat in the application of the presented methods is that the form of the constitutive model is assumed to be true for the purposes of parameter estimation and uncertainty quantification. This assumption carries an important limitation: while the four-fiber-family model used herein is microstructurally motivated, it is nevertheless phenomenological, thus violating the above assumption by definition. (We note that no truly mechanistic models currently exist to describe vascular mechanics, so this limitation is unavoidable.) Nevertheless, using our current model, we have shown that tissue-level quantities of interest such as tangent stiffness can be estimated with fairly high credibility for a variety of specimen types even when the model parameters themselves have higher uncertainties and lack physical interpretation. Caution is warranted in cases of non-pathological (e.g., development) or pathological (e.g., thrombus, calcification) situations wherein the model may have a more limited ability to reproduce the data. Model misspecification introduces an additional layer of uncertainty and potential bias to the analysis, separate from the quantified uncertainties in the parameter values. 46 It is therefore crucial to confirm that the adopted model has adequately captured the observed mechanical response-for example, by requiring a high coefficient of determination in regions that are included in the analysis, 5 and by qualitatively verifying independence of the residuals. As an extension of the present work, Bayesian model selection methods can be employed to choose between competing constitutive model forms. 23,47

CONCLUSION AND FUTURE WORK
In summary, we integrated a novel uncertainty quantification and propagation pipeline within our inverse modeling approach to estimate local mechanical properties in vascular tissues. The methods rely primarily on well-established empirical and analytic Bayesian techniques, and facilitate both subject-specific and group-level statistical analyses of the mechanical response of the vessel wall. Future directions include the application of the presented pipeline to mechanobiological network-based models of cell signaling as well as predictive models of vascular development, growth, and remodeling. [29][30][31] Other potential applications in disease include the comparison of aneurysmal, tortuous, dissected, and/or thrombotic regions of the vessel wall to their normal counterparts. 7,48 The approach can also be extended by accounting for uncertainties in the acquired data, including histologically prescribed tissue constituent mass fractions, OCT-derived local thickness estimates, and experimental vessel deformations relative to the unloaded configuration. 3,49 FIGURE A1 Automatic identification of fiducial markers in panoramic images. An original acquired image is binarized using an intensity threshold determined via Otsu's method, 18  proper calibration of the imaging system and subsequent mapping of the markers' image coordinates to real-world positions.
Intensity thresholding of the acquired image via Otsu's method 18 yields a binarized image in which the fiducial markers as well as most of the tissue region, the corners of the field of view, and spurious dark pixels display as black ( Figure A1). At this stage, segmentation of the binarized image into marker versus non-marker regions previously required interactive input through a graphical user interface. We automated the marker identification process by first computing the overall area of each dark region as well as its circularity, which is proportional to the ratio between its area and the square of its perimeter. Since the number of markers greatly exceeded the number of non-marker dark regions, the dark region with the median area and circularity could be classified as a marker with near certainty. Moreover, because all markers had a similar size and shape, all the remaining markers could be expected to have area and circularity close to the median values, relative to the overall range of observed area and circularity values. Each dark region was therefore categorized as a marker if its [circularity, area] vector was one of the 245 nearest neighbors to the median [circularity, area], as measured by the Mahalanobis distance, 50 which is computed using the standardized (i.e., normalized by standard deviation) circularity and area values. To establish correspondence between images from different views, the markers were then automatically sorted based on their azimuthal angle and distance from the mean of the marker centroid locations.
The correlation-based local image registration methods which are the foundation of pDIC depend on a high degree of similarity in the horizontal and vertical directions between the two images of interest; the acquired panoramic images, however, which correspond to eight radially symmetric camera positions, exhibit circumferential and radial similarity instead. Therefore, as in previous works by Genovese et al. 4,9 and Bersi et al., 5,6 pDIC was performed not using the acquired images directly, but instead on "rectified"/"unwrapped" images; that is, an interpolation of the original image data on a polar coordinate system. Because the camera is not exactly coaxial with the vessel and conical mirror being imaged, the center point of the acquired image does not correspond in general with the vertex of the conical mirror about which the image data should be unwrapped. Unwrapping the image data about a point far from the conical vertex yields characteristic undulation patterns in the interpolated image data ( Figure A2), which can cause pDIC to perform poorly. Previously, the search for an acceptable unwrapping center point relied on several iterations of interactive guessing, until an unwrapped image with minimal distortions was obtained; this process was required to be performed for every image acquired (336 images total per vessel, if the examined configurations correspond to the protocols described in Section 2.2.2).
For the present study, we automated this procedure using the following approach: First, the image data is unwrapped using the center of the image as an initial guess of the true center of the imaging setup. Under an optimal unwrapping, the inner boundary of the fiducial marker strip, which is nearly circular in the original image, would map to a nearly straight and horizontal line in polar space. The location of the marker strip boundary can be detected column-wise in the unwrapped image using an intensity threshold, since the marker strip is much brighter than the tissue region ( Figure A2). Due to the inaccuracy of the initial guess, the resulting unwrapped marker strip exhibits the undulation pattern described above, where the vertical coordinate of the strip boundary corresponds to the radius coordinate of a circle not centered at the origin. Specifically, where is the radius coordinate in the (sub-optimally) unwrapped image, is the angular coordinate, 0 is the true constant inner radius of the marker strip, and [ , ] is the discrepancy between the true center point and the initial guess. Using this relation, the optimal center point can be identified in one iteration by performing a nonlinear least-squares regression of the marker strip radius curve with respect to 0 , and ( Figure A2). The original image is then unwrapped with respect to the optimal center point, so that it may be precisely registered to corresponding images from different views and tissue configurations using pDIC.
Using methods described in detail by Genovese et al., 4,9 we used the pDIC-based image registration results to solve for the real-world 3D position of each point in the acquired images, which yielded a point cloud that represented the vessel wall geometry with high fidelity ( Figure A3a). Computation of the local deformation gradient field requires a surface representation of the vessel wall, however, which must be smooth to allow for differentiation of the deformed material point coordinates with respect to their coordinates in the reference configuration. Previously, the point cloud representation was converted to a surface through the use of regression using non-uniform rational basis splines, with the final result depending on subjectively selected shape parameters to compromise between surface smoothness and deviation from the point cloud data. To remove the need for user input, we developed the following automated procedure, which yields a fully parametrized smooth surface reconstruction of the vessel wall in the reference configuration (in vivo axial stretch under 70 mmHg pressure): First, the vessel centerline is approximated by applying a Gaussian filter to the and point cloud data in the direction (i.e., along the central axis of the vessel), with the

FIGURE A2
Automatic polar parametrization of panoramic images. For each panoramic image, the image data is interpolated on a polar coordinate system ("unwrapped") using the center of the image as an initial guess of the true center of the imaging setup. Inaccuracies in the center point about which the image is unwrapped will cause obvious undulations in the array of fiducial markers used to map the image data to real-world coordinates. At each angular coordinate value in the polar image, the radius of the boundary of the fiducial marker strip can be identified using an intensity threshold. Since the strip boundary ( ) will closely match the radius of a translated circle (i.e., one whose center is not the origin), an optimal correction vector [ , ] to the initially guessed center point can be determined via least-squares regression. The original image is then unwrapped with respect to the optimal center point, so that it may be precisely registered to corresponding images from different views and tissue configurations filter span equal to 5% of the total vessel length. Each point is then parametrized with respect to a polar cylindrical coordinate system centered locally at the corresponding centerline location ( Figure A3a). To arrive at a smooth surface representation, the local radius coordinate field is modeled as a Gaussian process, 51 with and as predictor variables ( Figure A3b). Regression is performed on data that are repeated in the direction, to enforce periodicity. The reconstruction and parametrization is completed by mapping the Gaussian process representation of the local radius coordinate field back to Cartesian space ( Figure A3c).
To transform the surface reconstruction of the reference configuration to all other configurations, the pointwise , , and displacement fields (computed using the point cloud data) were similarly modeled as Gaussian processes, and the resultant smooth displacement fields were applied to the reference geometry.

FIGURE A3
Smooth surface reconstruction of the entire vessel wall. (a) Starting from a panoramic digital image correlation (pDIC)-derived point cloud of the vessel wall, the vessel centerline is approximated by applying a Gaussian filter in the direction to the and data. Each point is then parametrized with respect to a polar cylindrical coordinate system centered locally at the corresponding centerline location. (b) To arrive at a smooth surface representation, the local radius coordinate field is modeled as a Gaussian process, with and as predictor variables. Regression is performed on data that are repeated in the direction, to enforce periodicity. (c) The reconstruction and parametrization is completed by mapping the Gaussian process representation back to Cartesian space