Feasibility of Vascular Parameter Estimation for Assessing Hypertensive Pregnancy Disorders

Hypertensive pregnancy disorders (HPDs) increase health risk for mothers and infants. It is estimated that 5% of pregnancies in the United States are complicated by pre-eclampsia, the most common form of HPD, and rates have been steadily increasing per year [1]. However, the pathophysiology of HPD is poorly understood, leaving clinicians without a screening test to reliably assess risk of adverse pregnancy outcomes [2]. While several studies have investigated the utility of biochemical markers and/or ultrasound (U.S.) in early prediction of HPD, their positive predictive value has not been high enough for standard clinical practice [3]. Even if HPD could be accurately predicted, there are still limited treatment options for HPD due to underinvestigation of the nature of the disorder [4].

It has been suggested that insufficient spiral artery remodeling during placental development may be the cause of high blood pressure and other hemodynamic disturbances in HPD [5]. Since the latest noninvasive in vivo imaging technologies have not been able to resolve the small spiral arteries (diameters of the order of ∼10² μm, which is currently smaller than the resolution of magnetic resonance imaging (MRI) for body imaging), a popular alternative research tool has been the Doppler U.S. velocimetry of the uterine arteries (UtAs), which are upstream from the spiral arteries and easier to visualize because of their larger diameter [6]. More recently, 4D flow MRI velocimetry of the UtAs has been investigated [7,8]. Compared to U.S., 4D flow MRI has lower temporal resolution and longer scan time but larger spatial coverage, improved spatial resolution, and measures velocity in three dimensions rather than only one.

In normal pregnancy, the UtAs remodel outwardly, with an increase in vessel lumen area and little to no change in wall thickness [9,10]. This is mediated by a combination of local wall shear stress, nitric oxide release, and local/systemic endocrine signaling. Animal studies on tissue reorganization within the UtA wall have observed hyperplasia and hypertrophy of the smooth muscle cells to maintain wall thickness during vasodilation [10,11]. Changes in elastin and collagen content in the UtA wall have been found to vary in different animal species [10]. One study reported evidence of increased myogenic tone in human pregnancy, which is surprising given the amount of vasodilator molecules released [12]. Nevertheless, this appears to signify that the uteroplacental arteries play a critical role in regulating the distribution of blood flow during pregnancy.

Doppler U.S. velocimetry has shown that UtA flow increases with gestational age (GA) while pulsatility and resistivity indices (PI and RI) decrease with gestational age in normal pregnancy [6,13]. The PI corresponds to the difference between the maximum and minimum blood velocity, normalized by the mean velocity and the RI corresponds to the difference between the maximum and the minimum blood velocity normalized by the maximum blood velocity. In early pregnancy, the velocity waveform may contain a diastolic notch, representing the reflected blood flow of uteroplacental circulation [14], which is expected to diminish by the third trimester as the uteroplacental vasculature continues to remodel [13]. These features suggest that persistently high pulsatility and the presence of the diastolic notch would be associated with adverse pregnancy outcomes. However, despite numerous investigations, U.S. measurement of UtA velocity still has low positive predictive value, which has precluded it from entering routine clinical practice [3].

Computational fluid dynamics has been a useful research tool to help identify sources of UtA velocity waveform changes and perhaps shed light on any biases that may weaken its diagnostic and screening utility. The earliest models were zero-dimensional (0D) lumped parameter models based on Frank's Windkessel model, which allowed mathematical descriptions of hemodynamic parameters using physical principles adopted from electrical circuits [15]. Fundamentally, the Windkessel model relates pressure (P) to flow (Q) as P = QR, where R denotes resistance. One-dimensional (1D) models have also been found to be computationally efficient and validated for accuracy in various arterial networks and flow conditions [16,17]. More recently, advances in computing have enabled three-dimensional fluid–structure simulations of blood flow in localized organs and tissues [18–20].

Simulations of UtA flow have reported that high resistance and small UtA radius recapitulated the abnormal waveform shapes (high PI with diastolic notch) that closely match measurements acquired from a pre-eclamptic mother [21,22]. Computational fluid dynamics methods were also leveraged to investigate spiral artery remodeling as a potential source of the abnormal waveform shapes detected in UtA U.S.. While insights were gained in the interaction between trophoblast invasion and hemodynamic changes of spiral arteries during pregnancy [5,23], evidence suggests that PI and the diastolic notch may be a moderate proxy for trophoblast invasion but other factors such as radial artery and arteriovenous anastomosis remodeling may also contribute to the waveform shape [24]. This suggests the need for more investigation into the complex remodeling-hemodynamic interactions underlying HPD that would lead to more effective clinical biomarkers than PI/RI. Previous uterine artery computational modeling has sought to understand the physical mechanisms underpinning measured ultrasound velocity wave-forms [25], but without considering a patient specific large vessel structure in the female pelvis.

In this study, we leverage a hybrid 0D–1D model first proposed by Sherwin et al. [26] to derive parameters, such as arterial resistance, compliance, and cross-sectional area, from noninvasive imaging flow measurements of the proximal uterine arteries for patient specific pelvic geometries. To achieve this, we developed a computational framework that combines a reduced-order Navier–Stokes model as a forward evaluation model and a Markov Chain Monte Carlo (MCMC) algorithm for predicting these biological indices that could be later used to assess the progression of HPD. We first demonstrate feasibility in the aorta with synthetic and measured MRI flow data, followed by a similar analysis of the maternal pelvic arterial network of a healthy pregnant subject and a pregnant subject that went on to develop hypertension later in pregnancy after MRI. We also compare the corresponding PI and RI from MRI and U.S. in these subjects.

In a Bayesian estimation setting, reduced-order models are commonly coupled with MCMC algorithms for inferring posterior distributions over the unknown parameters given noisy measurements. For example, Larson et al. [27] introduced a Bayesian model selection algorithm for identifying the position of wall abnormalities in a complex arterial network. Their algorithm employed a transition MCMC strategy coupled with a reduced-order flow model in order to examine different flow model setups, chose the one that best described the underlining phenomena, and matched the target waveform. Thus, the objective of this work was to solve a classification problem, where the boundary conditions of the problem were assumed and the geometrical and structural properties of the vessel were varying. Moreover, Colebank et al. [28] proposed a strategy that employed a delayed rejection adaptive metropolis algorithm in order to infer the boundary conditions for a three-element Windkessel model in pulmonary hypertension and performed a sensitivity analysis to identify the model parameter influence on the predicted pressure. In this case, a parameter reduction strategy was employed and nominal parameter values for R and C were assumed, which then reformed the inverse problem to identifying scaling multipliers for the nominal parameters at each outlet. The authors employed pressure data for the main pulmonary artery in order to infer the arterial resistance and compliance for the outflow vessels and performed model sensitivity analysis to assess the parameters' importance to the model output. In another approach, Păun et al. [29] proposed an MCMC method for inferring flow parameters in a pulmonary circulation vascular network in mice. The strategy was based on delayed rejection adaptive metropolis, with an additional parameter scaling technique. Moreover, they proposed starting the algorithm using a Maximum Likelihood approach in order to speed up the convergence of the algorithm and presented convergence tests in order to assess the effectiveness of the algorithm. In this work, we perform direct inference of the model parameters only via sparse, cross-sectional averages of 4D flow MRI measurements. This increases the complexity of the inference, as there is an inherent uncertainty caused by the discrepancy between the chosen model and the real flow conditions. Because of this discrepancy, different combinations of inputs might provide results that are in close agreement and thus provide high values of the likelihood function during the Monte Carlo sampling.

In summary, the main contribution of this work is a patient-specific methodology for directly inferring the values of arterial resistance, compliance, and equilibrium cross-sectional area, as well as an estimation of the absolute pressure in maternal pelvic arteries. This is a first step toward investigating the relationship between vascular parameters and pregnancy outcomes for early prediction of HPD.

For this study, we use the same aorta dataset reported in Kissas et al. [30]. The MRI data were acquired from a 27-year-old healthy female volunteer using a 1.5 T MRI scanner (Avanto, Siemens Healthineers, Erlangen, Germany). The protocol included a balanced steady-state free precession (bSSFP) localizer, prospectively electrocardiogram-gated two-dimensional (2D) cine, and 2D phase contrast images prescribed at four locations in the aorta and one location in the left common carotid artery. The temporal resolution was 29.4–34.8 ms and 20.7 ms for the 2D cine and 2D phase contrast MRI, respectively.

Using the same scanner, MRI data were acquired from a healthy pregnant volunteer, denoted as the normal subject (maternal age = 24 years, gestational age (GA) = 18.7 weeks) and a pregnant volunteer that developed hypertension in late pregnancy, denoted as the pre-HPD subject (maternal age = 25 years, GA = 19.1 weeks). The pre-HPD subject was neither diagnosed with chronic hypertension nor was on any related medication and was diagnosed with pre-eclampsia four months after the MRI data were collected. Each subject was positioned feet-first supine with a 12-channel spine array coil and two four-channel body matrix coils. Both patients were imaged in supine position with no tilt. An electrocardiogram monitor was attached to each subject for synchronization of the MRI to heart rate. The protocol consisted of a half-Fourier acquisition stimulated echo localizer of the abdomen and pelvis, a time-of-flight (TOF) angiogram of the abdomen and pelvis, a 2D prospectively gated phase contrast scan of the descending aorta, and prospectively gated 4D flow MRI of the uterine arteries (UtAs) and external iliac arteries (ExtIls). The details of the MRI acquisition parameters are listed in Table 1.

During postprocessing of the MRI data, the TOF angiograms were segmented (Seg3D [31]) and centerlines were extracted (VMTK [32]) to compute path length from the descending aorta to the UtAs and ExtIls (Fig. 1). The equilibrium areas of each vascular segment were estimated from multiplanar re-formats of the TOF images. Velocity and area wave-forms were extracted from 2D phase contrast images using ImageJ. The 4D flow images were processed with custom software (MATLAB, Natick, MA) [33], followed by velocity-based thresholding of the volumetric isosurface (Ensight, CEI; Apex, NC). Velocity wave-forms were extracted from the left/right UtAs and left/right ExtIls. Since both the 2D phase contrast MRI and 4D flow MRI were prospectively gated, a few end-diastolic cardiac phases were not acquired from each cardiac cycle. These end diastolic velocity values were extrapolated by averaging the velocities of the first and last cardiac phase (Table 1).

Each subject was positioned in semirecumbent supine during U.S. scanning by a clinician (N.S.) with extensive experience in prenatal U.S.. The UtAs were scanned using the C4-8 transabdominal probe of a GE Voluson E10 (GE Healthcare, Wisconsin) U.S. machine. Measurements of UtA Doppler velocity wave-forms were recorded bilaterally, and PI and RI were calculated for comparison with the MCMC estimation. Pregnant subjects were enrolled in this study approved by the institutional HIPAA-compliant Institutional Review Board at the University of Pennsylvania.

where $p(x,t)$⁠, $q(x,t),$and $A(x,t)$represents the pressure, the flow and cross-sectional area, respectively. Moreover, $β=π h0E(1−ν2) A0$ and $A0$ denotes the vessel's cross-sectional area at equilibrium, E the Young's modulus, h₀ the wall thickness, ν the Poisson ratio, and ρ the blood density. Values for E and ν can be found in the literature, while $h0$ and $A0$ can be measured via MRI. Moreover, $KR$ is a friction parameter which depended on the chosen velocity profile and α is a momentum flux correction parameter, which accounts for the nonlinearity of the sectional integration in terms of the local velocity [26]. In this study, $KR=−22 μπ$⁠, where $μ$ is the dynamic blood viscosity and $α=1.1$⁠. For all experiments, $β$ is computed based on the empirical relation, between the Young's modulus, the equilibrium cross-sectional area and the vessel wall-thickness reported by Olufsen et al. [35]. The tortuosity of the geometry is assumed to be small enough such that the system (1) is valid, but in reality there was some uncertainty induced by neglecting the complex topology, particularly of the uterine arteries.

where $p1,p2,3, (ûx)1, and (ûx̂)2,3$ are the pressure and velocity values of the parent and daughter vessels, respectively.

This constitutes the full forward model for predicting the target quantities (flow, wall displacement, and absolute pressure). We employ an in-house discontinuous Galerkin (DG) solver [26,37] to compute the solution of this model. We extract the vessel centerlines from the three-dimensional arterial geometry of each patient and created a graph as the input geometry of the DG simulator.

for each vessel j.

where N_o is the number of outlets, $q̂j$ represents the target measurements, $qj$ represents the forward model prediction for parameter set $ϑj$⁠, and $τj$ represents the noise variance of the measurements for vessel #j. The use of the Gaussian likelihood is based on the assumption that the measurements may be corrupted by Gaussian noise with zero mean and variance $τj−2$⁠.

where $ζj, ηj$ are user-defined lower and upper bounds, respectively, of the parameter range of the considered outlet. Furthermore, $γj$is the lower and $δj$the upper bounds for the noise parameter for each outlet. The uniform distribution is chosen for the priors [40] to limit the chance of having the distribution mass concentrated to a particular region. In the case that information over the prior distribution form and range is known, informative priors can be employed in order to accelerate the sampling process convergence.

where N denotes the number of posterior samples. The above algorithm provides a statistical description of the predictions, which later is used to assess the quality of our conclusions.

Flow rates are computed by multiplying the predicted velocity wave-forms with the predicted cross-sectional area for each vessel. Then, they are compared to the target flow rates, which are computed by multiplying the cross-sectional area measured from TOF MRI with the velocity measured from 4D flow MRI. The algorithm is schematically presented in Fig. 2.

The Bayesian inference algorithm was tested on three MRI volunteer datasets, one collected in the aorta and two in the maternal pelvic arteries (Sec. 2.1). For all experiments, we predefined parameters $KR=−22μπ$⁠, $α=1.1$ to account for viscous losses, the blood density $ρ=1060 kg/m3$⁠, and the viscosity $μ=3.5 mPa s$⁠. We provide a sketch of the vascular networks considered in this paper, as well as the connection between the different vessels and the boundary conditions in Fig. 3.

The code for all the methods is created by combining the python library for Bayesian statistical modeling PyMC3 [41] with a custom C language code including a custom loss function in the python library Theano [42].

The arterial network geometry consists of four vessel segments with one inlet (ascending aorta), a bifurcation, and two outlets (left common carotid artery (LCC) and descending aorta (dAo)), as shown in Fig. 1. Given this dataset, we perform two computational experiments, one using synthetic velocity data and one using MRI measured velocity data. This aorta dataset was originally reported in a previous study, see Ref. [30].

In the first aorta experiment, the geometry was coupled with synthetic velocity data at the outlets to formulate a simplified problem as a proof-of-concept of the proposed Bayesian inference framework. To construct the synthetic data, the DG simulator was provided with a set of Windkessel parameters (R, C). The generated velocity data was then corrupted with 2% Gaussian noise to simulate MRI measurement error. Despite neither the noise level being that small nor the noise model Gaussian in real cases, we consider this case as a sanity check, thus we define it in a simple manner. Moreover, we consider the model discrepancy [43], the difference in approximation quality between the employed and the real physical models that describe the flow, as the main factor contributing to the prediction uncertainty. The period of the volunteer's cardiac cycle T = 0.78 s was preserved in the simulation and $Np=50$⁠, where $Np$ the number of measurements, were chosen to represent the waveform. Further details of the simulation setup have been previously discussed in Ref. [30].

The Bayesian inference algorithm is setup to take the geometry and synthetic data and estimate the Windkessel parameters and pressure wave-forms. The search space is defined by uniform prior distributions for the unknown parameters (see Table 2) chosen such that we do not assume any prior knowledge on the distribution mass, as explained in Sec. 2. The area wave-forms provided by MRI are not taken into consideration for simplicity of this proof-of-concept experiment, but MRI measurements of equilibrium cross-sectional area are used to determine the wall stiffness parameter β in the one-dimensional blood flow model. To assess accuracy after convergence, the estimated parameters are compared with the previously validated [30] parameters $RdAo=1.667×108 Pa s m−3, CdAo=9.002×10−9 Pa−1 m3, RLCC=2.102×109 Pa s m−3$⁠, and $CLCC=2.538×10−10 Pa−1 m3$⁠, which are taken based on the findings reported in Ref. [30].

In the second aorta experiment, the same geometry is coupled with the measured MRI velocity data at the outlets in order to estimate the set of Windkessel parameters that provide the best match for the given previously validated measurements. Prior to starting the Bayesian inference algorithm, a Gaussian process regression [36] with a periodic kernel is performed on the MRI wave-forms to ensure periodicity, to augment the $Np=38$ points to $Np=100$ to improve precision and to align the predicted velocity and MRI velocity wave-forms which have different temporal resolutions and starting times. More details about the preprocessing of the MRI data are described in Ref. [30]. The Bayesian inference algorithm is setup with a search space defined by the uniform prior distributions reported in Table 2.

The Bayesian inference methodology is performed on two data-sets comprised of MRI measurements of velocity from the female pelvic region of two human volunteers. The first data-set corresponds to the normal subject and second to the pre-HPD subject, see Sec. 2. The maternal pelvic arterial network geometry consists of seven vessels extending from the descending aorta (dAo) to the right and left common iliac arteries, to the right and left uterine arteries (RUtA, LUtA) and external iliac arteries (RExtIl, LExtIl). Each network contained one inlet, four outlets, and two bifurcations, as shown in Fig. 1(c) and 1(e) for the normotensive and the pre-HPD subjects, respectively. For both cases of maternal pelvic geometries, vessel 3 corresponds to the Aorta, vessel 6 to the Right Common Iliac artery, vessel 2 to the Left Common Iliac artery, vessel 7 to the Right Uterine artery, vessel 4 to the Left Uterine artery, vessel 1 to the Left External Iliac artery, and vessel 5 to the Right External Iliac artery.

For each network, the Bayesian inference algorithm is provided with MRI measured velocities at each outlet and then executed to estimate the Windkessel parameters R and C as well as the equilibrium cross-sectional area (A₀), similar to the aorta experiments. Unlike the aorta experiments, cross-sectional area measurements are included in the inference. In contrast to the descending aorta region, where the cross-sectional area magnitude of the vessels is large enough to be accurately captured by the MRI resolution, in the pelvic region the arteries are very narrow and thus inaccuracies are introduced in the measurements. In the maternal pelvic experiments, the structural MRI measured cross-sectional area of the vessels that do not correspond to outlets (i.e., descending aorta, left and right common iliac arteries) is inserted to the model and kept constant throughout the entire inference procedure, to avoid identifiability issues. Therefore, we consider the cross-sectional area as a free parameter for the outlets to potentially correct for the induced inaccuracies in the geometry of their parent vessels.

Since 4D flow MRI consists of velocity measurements from planes along the vessel cross section at different locations, we can construct time series for both the maximum and the mean velocity. In standard practice, the mean velocity wave-forms are used to compute the mean flowrate. However, the mean velocity wave-forms have low pulsatility therefore the maximum velocity is used in practice for calculating PI and resistivity index (RI). The reduced-order model employed in this work is a cross-sectionally averaged version of the Navier–Stokes equations, so we will use the mean MRI velocity for estimating the pressure wave-forms, R, C, and A₀. Moreover, we will run a separate Bayesian inference experiment for the sole purpose of calculating PI and RI to compare with U.S. measurements. This computation constitutes another sanity check for assessing the quality of the approximation that the framework can achieve for the case of predicting biomarkers. The biomarkers are computed using the mean, the minimum, and the maximum of the velocity waveform; therefore, we can assess how well the model can approximate these waveform characteristics.

For both volunteers, the measured 4D flow MRI velocity waveforms have a low temporal resolution, consisting of only $Np=14$ points over a cardiac cycle of T = 0.612 s. To this end, the MRI velocity wave-forms are augmented by a Gaussian process with a periodic kernel regression to create a set of $Np=55$ points. This additional step is beneficial to the procedure because by creating a larger baseline periodic dataset we can have stronger comparison with the prediction which results to a more accurate algorithm. Moreover, this step is used to align the velocity time-series produced by the model and the MRI measurements so that the predictions and the measurements are on the same times. Table 3 contains the bounds of the uniform prior distributions predefined for the Bayesian inference framework employed for the normal subject and the pre-HPD subject.

The pelvic arterial network parameter estimation simulation takes around 5 days on a single core of a Thinkpad p52 s laptop with an Intel Core i7-8650 U CPU and 32GiB of RAM. This computational cost could be naively cut down by considering a more powerful workstation and less steps of the sampling algorithm. We considered more steps than we could in order to be sure that the sampler has converged. For more patients, we can asynchronously parallelize the process, meaning that we run multiple processes for multiple patients at the same time without waiting for one process to finish to start the next one.

We present the velocity and the pressure wave-forms for the case of the synthetic data in Fig. 4 and the predicted parameters in Table 4. In the aorta model with synthetic velocity data, the predicted outlet velocity wave-forms closely matched the target synthetic wave-forms, see Figs. 4(a) and 4(b). We draw 200 samples of the posterior parameter distribution and calculated the flowrate using the predicted velocity and area. For this case, the flowrate resulting from the parameters with the maximum probability closely matched the target. For the descending aorta the model prediction is 4577.9 ± 8.1 and the target 4578.8 mL/min which results to −0.02% error. For the left common carotid artery, the predicted flow is equal to 462.8 ± 6.8 and the target 466.6 mL/min which results to −0.82% error. The predicted Windkessel parameters also closely match the target values for the left common carotid (resistance $R=2.08×109$ versus $2.10×109 Pa s m−3$⁠, −1.09% error; compliance C = 2.54 × $10−10$ versus 2.54 × $10−10 Pa−1 m3$⁠, 0.12% error) and descending aorta (R = 1.67 × 10⁸ versus 1.67 × 10⁸$Pa s m−3$⁠, 0.24% error; compliance C = 8.981 × $10−9$versus 9.00 $10−9 Pa−1 m3$⁠, −0.28% error) (Table 4). In addition, the aorta model provided absolute pressure wave-forms for both outlets (see Figs. 4(c) and 4(d)). Although validation waveforms are not available from the measured MRI data, the predicted pressure time-series are within physiological range. The absolute systolic/diastolic pressure is 115.2 ± 1.4/59.7 ± 1.5 mmHg in the descending aorta and 108.0 ± 1.4/61.8 ± 1.5 mmHg for the left common carotid artery.

When the MRI measured velocity data are used in place of synthetic data, the aorta model still achieved close agreement with the target MRI measurements. The predicted flowrate in the descending aorta is 4803.9 ± 17.5 versus 4519.1 mL/min, 6.3% error and in the left common carotid artery 353.3 ± 6.8 versus 343.7 mL/min, 2.8% error (Figs. 5(a) and 5(b)). Although ground truth the resistance and compliance values are not known for this particular dataset, the inferred parameters can be determined the descending aorta resistance is 1.45 × 10⁹ Pa·s m⁻³ and compliance is 1.88 × 10⁻¹⁰Pa⁻¹ m³. The predicted resistance and compliance in the left common carotid artery are 1.15 × 10⁸ Pa·s m⁻³ and 1.29 × 10⁻⁸ Pa⁻¹ m³, respectively, see Table 4. The predicted pressure wave-forms are also within physiological range [44]. The systolic/diastolic pressures in the descending aorta are 102.5 ± 1.0/55.1 ± 1.0 mmHg and in the left common carotid artery pressure are 95.4 ± 0.7/55.1 ± 0.8 mmHg (see Figs. 5(c) and 5(d)).

The computational model proposed in Sec. 2 for estimating the parameters in both the normal subject and the pre-HPD subject cases. The predicted flowrate in the left and right external iliac arteries across 200 samples of the posterior distribution in the normal subject, are 950.18 ± 14.04 and 919.64 ± 38.92 mL/min, respectively. The predicted flowrate in the left and right uterine arteries are 323.02 ± 3.84 and 65.57 ± 0.98 mL/min, respectively (Fig. 6). The pre-HPD subject has left and right external iliac artery flow rates 472.2 ± 16.5 and 711.02 ± 27.8 mL/min, respectively, and left and right uterine artery flowrate of 62.75 ± 4.5 and 134.22 ± 4.1 mL/min, respectively (Fig. 6). When comparing between predicted and MRI-measured flow rates, the errors are highly variable (5.29–92.4% error). These discrepancies occurred because the model cannot provide accurate predictions of the cross-sectional areas of the right and left external iliac arteries. One reason for this lack of predictability is high errors are accumulated for the case of the cross-sectional area, as it is present in the denominator of the β parameter in the pressure equation of the model. Another reason is that we assume that the area measurement provided by the structural MRI is the equilibrium cross-sectional area for the nonoutlet vessels, left and right common iliac arteries and the descending aorta, which can contribute to inaccuracies. We present the predictions for the biomarkers computed from MRI, Bayesian inference and, also, include the Doppler Ultrasound derived biomarkers on Table 5. The predicted PI and RI values of the Normal subject are collectively lower than that of the pre-HPD subject. This trend agrees with the corresponding MRI and ultrasound-based PI and RI comparisons.

The predicted R and C values are reported in Table 6. The predicted R values for the pre-HPD subject are higher than the R values of the Normal subject in the left external iliac artery, left uterine artery and right external iliac artery. The predicted compliance for bilateral external iliac arteries is lower in the pre-HPD subject compared to that of the normal subject. However the opposite is true for bilateral uterine arteries, where the compliance is higher in the pre-HPD subject. Table 6 also reports the estimated cross-sectional area values. Both subjects have similar external iliac artery areas with the pre-HPD subject having a slightly narrower left external iliac artery. For the pre-HPD subject, both uterine artery cross-sectional areas are close to 25% of the external iliac artery areas. For the normal subject, we see a higher variability in the uterine cross-sectional areas with the left uterine artery higher. In terms of the predicted pressure wave-forms for both cases, they are (in general) higher for the pre-HPD subject case for systolic and diastolic pressures compared to the Normal subject case. In the left uterine artery, the predicted pressure is 152.4 ± 7.8/65.1 ± 7.3 mmHg versus 107.78 ± 1.83/32.99 ± 1.49 mmHg, and in the right uterine artery, the predicted pressure is 148.2 ± 8.3/64.7 ± 7.14 mmHg versus 85.78 ± 3.33/31.0 ± 1.84 mmHg (see Figs. 6 and 7). Moreover, from the comparative figures for the velocities and the pressures in Figs. 6 and 7, respectively, we observe that the velocities are collectively higher in the case of the normal subject while the pressures are collectively higher in the case of the pre-HPD subject.

This study demonstrates the feasibility of inferring parameters of the maternal pelvic arterial network from MRI geometry and velocity data within a computational blood flow modeling framework. By adopting a Bayesian inference approach, we are able to iteratively match the predicted velocity wave-forms with the measured velocity wave-forms in a favorable manner. By doing so, we are able to recover parameters such as resistance, compliance, and pressure, which are not measured but inferred from physical principles. We initially perform proof-of-concept experiments in an aorta and carotid dataset of a healthy volunteer and show that the Bayesian inference algorithm matched the target synthetic flow wave-forms and the target MRI flow wave-forms. There is greater uncertainty in the predicted velocity wave-forms for the aorta model with MRI flow data potentially because of the measurement noise present in the dataset. In the synthetic data case, the target resistance and compliance are chosen to be the same as the ones computed by a different estimation technique proposed by Kissas et al. [30], when applied to the real MRI descending aorta case. We choose these values because we are already aware that the results are in physiological ranges. Both the target and the predicted parameters are similarly matched in both outlets. In both aorta models, the left common carotid artery had lower predicted resistance and higher compliance than the descending aorta. The predicted pressure wave-forms of the descending aorta and left common carotid artery are also within physiological range based on cardiac catheterization measurements in a human [44].

In the maternal pelvic arterial models, the Bayesian inference algorithm successfully recovered the parameters despite greater complexity in the network. Traditional measures of PI and RI show that the subject who went on to develop pre-eclampsia (pre-HPD subject) has elevated second trimester uterine artery pulsatility compared to the normal case (normal subject), a trend which is often interpreted as higher resistance. When comparing the predicted resistance values between the models for the Normal subject and the pre-HPD subject, the pre-HPD subject has higher resistance than the normal subject in the left uterine, the right and the left external iliac arteries, which indicates potential incomplete remodeling of these vessels in this case. On the other hand, the right uterine artery has higher resistance in the normal subject potentially because it has not remodeled yet. Previous studies have described reduced vascular tone as part of the general pattern of arterial remodeling in pregnancy [10]. Since pre-eclampsia has been described to involve insufficient conversion of the spiral arteries from prepregnant vasoactive vessels to flaccid conduits [5], we postulate that HPD results from failure to also completely lower the compliance of the uterine arteries. However, our results showed the pre-HPD subject case had higher compliance than the normal subject case for the uterine arteries. Although the discrepancy between our results and this theory can be attributed to patient variability and the pilot nature of this initial study, animal studies have indeed found higher myogenic tone in healthy pregnant sheep and guinea pigs compared to healthy nonpregnant animals [10], which seems to show that there may be a complex relationship between resistance and compliance and clinical outcomes that needs to be further investigated. Moreover, the uterine artery is potentially exposed to more local vascular mediators from the placenta/pregnancy, possibly explaining why we may see remodeling/compliance changes there that are different from what is seen in systemic circulation [45].

Unlike the aorta, the outlet cross-sectional areas for the maternal pelvic arterial network were also considered as inferred parameters because the 4D flow MRI did not have sufficient spatial resolution to precisely measure the areas of the small uterine arteries. The pre-HPD subject has lower cross-sectional area in the left uterine artery compared to the normal subject, while higher in the right uterine artery. Given the hypothesis that smaller luminal areas correlated with poorly remodeled uterine arteries in pre-eclampsia [10,21], it made sense that the pre-HPD subject would have low cross-sectional areas of the uterine arteries, whereas the Normal subject would have at least one dilated uterine artery. For the right uterine artery, the predicted cross-sectional area of the normal subject is smaller than the pre-HPD subject, which might be an indication that the particular artery has yet to start remodeling. The pre-HPD subject had consistently higher predicted pressure in all four outlets compared to the Normal subject.

There are several potential reasons for the discrepancies observed between the predicted and target velocity wave-forms. First, the measured velocities from MRI have a significant amount of noise, so they do not necessarily obey the underlying physics of the Navier–Stokes equations. Second, we performed simulations on a single subject using a single measurement and therefore did not have the advantage of averaging velocities over a population which would have improved robustness. Third, we assumed that the curvature of the vasculature is small enough to be considered one-dimensional. We also, assumed that the velocity magnitude in the direction of the flow is much larger than in any other directions. These assumptions were made in order to reduce the model from three-dimensional to one-dimensional. Due to the complexity of real flows, these assumptions do not hold completely, therefore we have lost accuracy and introduced model discrepancy in the process. Finally, the low spatiotemporal 4D flow MRI resolution in the region of the maternal pelvic arterial network introduced uncertainties due to the inaccuracies in structural parameter measurements, such as the vessel cross-sectional area, and the small number of velocity measurements comprising a cardiac cycle. These modeling assumptions and measurement errors may result in phenomena where the uncertainty provided for the model might be high for some parameter of interest, i.e., see the right external iliac artery prediction in Fig. 6. These shortcomings of the method, which are related with the imaging techniques and technology, could be potentially alleviated by the advances in medical imaging technologies, paving the way for more accurate parameter estimation and potentially the successful enhancement of the proposed methodology in higher dimensions (i.e., three-dimensional flow data). Parameter inference with higher dimensional data is possible but not efficient at this point, one reason being that the lack of measurement accuracy lessens the convergence speed of the inference methods.

Granted that the goal of this method is to be applied in the clinic, we aim to speed up the Bayesian inference by addressing the wall-clock time bottleneck in our process. This process can be sped up as explained in Sec. 2.3 by considering some educated initial guess for the parameters. In this work, we considered no prior belief about the values of the parameters, meaning that we defined a huge space to search for the parameters that would provide us with velocity waveforms that closely match the 4D MRI results and no bias on where the algorithm should search. This means that it takes many iterations for the sampler to find the area where with high probability the resulting velocity will match the measured one. If we considered prior knowledge the algorithm would have probably converged faster.

Recently, a growth and remodeling framework was proposed by Gleason and Sedaghati [46] for the maternal pelvic arterial network of a normal pregnant subject. For the UtAs, we observed that the values of the discovered diameters (LUtA = 5.3 mm, RUtA = 2.67 mm) are within the bounds found in Fig. 8 of Gleason et al., as well as the discovered resistance values (LUtA = 0.6 × $105$ dynes s/cm⁵, RUtA = 0.1 × $105$ dynes s/cm⁵). The mean flow value for LUtA is 323 mL/min, which is within the values presented in Fig. 8 of Gleason et al., while the mean flow value for RUtA is 65 mL/min, which is less than the minimum value presented in the same figure. One potential source of this discrepancy is that the values of mean velocity wave-form for the Uterine arteries in Gleason et al., Fig. 10, is higher than the one measured from MRI in our paper. For the External Iliacs, we observed that the values of the discovered diameters (LExtIl = 9.3 mm, RExtIl = 9.3 mm) are larger than the upper bound found in Fig. 13 of Gleason et al. The discovered Resistance values (LExtIl = 0.4 × 10⁴ dynes s/cm⁵, RExtIl = 0.4 × 10⁴ dynes s/cm⁵) and the mean flow rates (LExtIl = 950 mL/min, and for RExtIl = 919 mL/min), are also higher than the values presented in Ref. [46]. A possible reason for the difference in the resistance values is that in our work we consider no relation between the two vascular parameters. Other possible reasons for the discrepancies between the two models are the underlining assumptions and the constitutive/phenomenological models considered for the different methods.

A limitation of this study was that the geometry of the maternal pelvic arterial models did not include all the branches of the internal iliac arteries aside from the uterine arteries. While it is known that a key feature of maternal physiological adaptation to pregnancy is diverting more blood flow to the uterine arteries [6], this simplification may have caused inaccuracies without including the other branches. A more comprehensive model would require improvements in the spatial resolution of the 4D flow MRI technique to measure velocities in the other branches or at least inclusion of their structures so that the additional outlet velocities can be inferred. Another limitation is that this study reported total peripheral resistance, rather than separating between characteristic impedance and the systemic peripheral resistance. This separation can be useful in future studies to determine the source of the resistance, whether in the larger arteries of the pelvis or the small uteroplacental vessels (e.g., arcuate, radial, spiral arteries), which can affect interpretation of the remodeling changes based on computational results. It should also be noted that two large vessels, the brachiocephalic trunk and left subclavian artery, are missing from the aorta geometry, which might affect the results of this study. However, the fact that our estimated pressures were within physiologic range is reassuring regarding the validity and applicability of our chosen methods.

This study demonstrated the use of Bayesian inference to predict parameters that are not routinely measured clinically but can be helpful in identifying high risk pregnancies in early gestation with the aid of MR imaging and computational fluid dynamics. In future work, more patients are needed to investigate the relationship between pressure, resistance, compliance, cross-sectional area, and disease characteristics. This will aid in determining the most influential parameters that can be correlated with clinical outcomes. These parameters are promising biomarkers because they may be more physiologically relevant compared to PI and RI.

We are grateful to Brianna Moon and Veronica Aramendia Vidaurreta for their help in acquiring the aorta MRI data. The maternal pelvic arterial network MRI data were acquired with the help of MRI technologists, research coordinators, and pregnant research participants at the Hospital of University of Pennsylvania. G.K. and P.P. would like to acknowledge support from the U.S. Department of Energy under the Advanced Scientific Computing Research program (Grant No. DE-SC0019116) and Air Force Office of Scientific Research (Grant No. FA9550-20-1-0060). E.H. would like to acknowledge support from the National Institute of Health (Grant No. 1F31HD100171). E.H., N.S., W.W., and J.D. acknowledge support from National Institute of Health (Grant No. U01HD087180). E.W.T acknowledges support from NIH Medical Scientist Training Program T32 GM07170.

dAo =

descending aorta

DG =

discontinuous Galerkin

HPD =

hypertensive pregnancy disorders

MCMC =

Markov Chain Monte Carlo

(L/R)CC =

(left/right) common carotid artery

(L/R) CI =

common iliac artery

(L/R)ExtIl =

(left/right) external iliac artery

(L/R)UtA =

(left/right) uterine artery

PI =

pulsatility index

RI =

resistivity index

TOF =

time of flight

U.S. =

ultrasound