J R Soc Interface

Computational modeling, Statistics, Applied mathematics

Mette S Olufsen

*Authors contributed equally.

# A computational framework for generating patient-specific vascular models
and assessing uncertainty from medical images

Michelle A Bartolo1,∗ Alyssa M Taylor-LaPole1,∗ Darsh Gandhi1,4 Alexandria
Johnson1,5 Yaqi Li1,6 Emma Slack1,7 Isaiah Stevens1 Zachary Turner1,8 Justin D
Weigand2 Charles Puelz2 Dirk Husmeier3 Mette S Olufsen1 1Department of
Mathematics, North Carolina State University, Raleigh, NC, USA
2Division of Cardiology, Department of Pediatrics, Baylor College of Medicine
and Texas Children’s Hospital, Houston, TX, USA
3School of Mathematics and Statistics, University of Glasgow, Glasgow, UK
4Department of Mathematics, University of Texas at Arlington, Arlington, TX,
USA
5Department of Mathematics and Statistics, University of South Florida, Tampa,
FL, USA
6North Carolina School of Science and Mathematics, Durham, NC, USA
7Department of Mathematics, Colorado State University, Fort Collins, CO, USA
8School of Mathematical and Statistical Sciences, Arizona State University,
Tempe, AZ, USA<EMAIL_ADDRESS>
###### Abstract

Patient-specific computational modeling is a popular, non-invasive method to
answer medical questions. Medical images are used to extract geometric domains
necessary to create these models, providing a predictive tool for clinicians.
However, in vivo imaging is subject to uncertainty, impacting vessel
dimensions essential to the mathematical modeling process. While there are
numerous programs available to provide information about vessel length, radii,
and position, there is currently no exact way to determine and calibrate these
features. This raises the question, if we are building patient-specific models
based on uncertain measurements, how accurate are the geometries we extract
and how can we best represent a patient’s vasculature? In this study, we
develop a novel framework to determine vessel dimensions using change points.
We explore the impact of uncertainty in the network extraction process on
hemodynamics by varying vessel dimensions and segmenting the same images
multiple times. Our analyses reveal that image segmentation, network size, and
minor changes in radius and length have significant impacts on pressure and
flow dynamics in rapidly branching structures and tapering vessels.
Accordingly, we conclude that it is critical to understand how uncertainty in
network geometry propagates to fluid dynamics, especially in clinical
applications.

###### keywords:

arterial network generation, change points, patient-specific modeling,
cardiovascular fluid dynamics, medical image analysis, uncertainty
quantification

††journal: rsif

##  1 Introduction

Medical images, such as computed tomography (CT) and magnetic resonance
imaging (MRI), are widely used for disease diagnostics, treatment planning,
and clinical studies [1, 2, 3, 4]. It is a non-invasive and reliable way to
assess the structure of the heart, lungs, and vasculature. Image segmentation
allows for increasingly precise analysis of anatomical data by isolating key
features of an image using computational software [5, 6]. However, there is
inherent uncertainty in the image segmentation process due to complex anatomy,
patient motion during scanning, variability in user expertise, and noise [7,
8, 9, 10]. Despite advances that include semi- and fully automatic extraction
using machine learning, most software requires manual editing to capture the
entirety of the desired anatomical features relevant for a given study [9]. In
particular, in vivo images are only captured up to a finite resolution and
certain anatomies, such as the rapidly branching structure of the pulmonary
vasculature, are nearly impossible to extract in an automated manner [11, 9].
Vascular networks span several orders of magnitude, but resolution within an
image is constant [12], adding significant uncertainty when extracting data,
especially for the smallest vessels [11]. Moreover, for images acquired
without contrast, differentiating pulmonary arteries, airways, and veins, and
identifying junctions is not a trivial task [9].

Combining medical imaging data with computational fluid dynamics (CFD)
modeling can provide insights into the structure and function of the
vasculature by predicting hemodynamic properties that cannot be measured in
vivo [13, 14]. Many studies use patient-specific geometry derived from medical
images to capture network dynamics in systems of branching blood vessels using
one-dimensional (1D) or three-dimensional (3D) CFD [15, 11, 16, 17, 18, 19,
20]. These studies superimpose dynamic pressure and flow data onto geometric
domains extracted from images and make predictions of blood pressure, flow,
shear stress, and wave intensity to simulate the impact of disease and
treatment. 3D models provide detailed predictions of local flow patterns in
regions with significant secondary flows such as in curved vessels or
junctions, but are computationally expensive [20, 21]. 1D pulse wave
propagation models make accurate and fast predictions, especially when there
are no sudden changes in cross-sectional area and regions of unidirectional
flow [22]. Several studies have compared 1D and 3D models and found reasonable
agreement in pressure and flow predictions at multiple locations throughout
the cerebral, coronary, and systemic vasculatures [23, 24, 25]. This study
utilizes a 1D arterial network model that propagates flow and pressure
predictions based on geometry from imaging data.

The process of network reconstruction and centerline extraction creates a 1D
characterization of blood vessel geometry from a 3D surface. Due to the
complex structure of vascular networks, there are multiple methods to generate
centerlines, potentially yielding different outcomes. Depending on centerline
definition, segmentation quality, and extent of smoothing applied, a given
vessel may be assigned different radii and length for each method, leading to
differences in hemodynamics. A popular centerline extraction program is the
Vascular Modeling Toolkit (VMTK) [26], which uses maximally inscribed spheres
to determine centerlines describing the length and radii of each vessel in a
network. Vessels widen as they bifurcate, causing spheres to intersect and
indicating a junction has occurred. This process is sensitive to intricate
geometries common in the vasculature, such as high tortuosity, fast branching,
and small vessels, which sometimes causes junctions to be placed outside of
physiological range. Ellwein et al. [27] note potential flaws in VMTK’s
centerline creation and use angles between vessels to adjust bifurcations.
SimVascular is another software used for centerline extraction, where users
create vessel centerlines manually [28]. This overcomes the issues associated
with maximally inscribed spheres, but is prone to user error and is time-
consuming for large networks [28]. Chen and Cohen [29] propose a method based
on geodesic or minimal paths to compute centerlines; however, this often
results in extracted paths being placed outside of the center of each vessel
and has only been tested in retinal networks. Skeletonization provides an
alternative option by iteratively removing voxels from segmented images until
a skeleton that preserves the original branching pattern is obtained [30, 31].
This is an effective way to create centerlines, but errors such as false
branches, double edges, duplicate points, and small cycles must be manually
corrected [32]. These studies demonstrate the need to investigate the network
creation process to ensure that vessels are assigned appropriate dimensions,
especially when used to gain medical insight [33].

Uncertainty arises in nearly every level of the modeling process, including in
the assumptions to derive the equations that model the physical system,
parameter choices, and geometry extracted from image segmentation. Uncertainty
quantification has been studied extensively in computational modeling of
hemodynamics. For example, Colebank et al. [11] study the influence of pre-
segmentation parameters, such as thresholding and smoothing, on 1D predictions
in ex vivo murine arteries. Chen et al. [34] incorporate intrinsic parametric
uncertainties to depict blood flow and pressure within a stochastic model of
arteries. Bertaglia et al. [35] highlight the effects of changes in elastic
and viscoelastic parameters in fluid-structure interaction models of arteries.
While these studies effectively quantify uncertainty in the modeling and
parameter estimation process, they do not examine uncertainties generated in
the network extraction process from in vivo human images. It is well known
that there is a significant interplay between vascular geometry and
hemodynamics, and therefore, uncertainty of geometric features of vessels must
be explored and quantified [36].

In this study, we investigate the process of extracting networks from medical
images and the associated hemodynamic uncertainty. Using 3D Slicer, we segment
the pulmonary arteries from a normotensive human CT scan and the aortic
vasculature of a double outlet right ventricle (DORV) patient from magnetic
resonance angiography (MRA). We generate multiple segmentations from each scan
– five for the pulmonary vasculature and ten for the aorta – to demonstrate
variations that occur between segmentations, despite using identical
thresholding parameters. In the pulmonary vasculature, each segmentation
captures a varying number of vessels, so to enable comparison, we standardize
networks using a radius pruning method. To overcome centerline extraction
uncertainty, we develop a novel algorithm to align networks and adjust
junctions to physiologically accurate locations within bifurcations, altering
the length and radius of each vessel. Subsequently, we employ change point
analysis to determine a representative radius for each vessel, avoiding
outlier data that may not belong to the given vessel (Figure 1). We propagate
uncertainty in segmentation, junction location, and radius choice to a 1D
fluid dynamics model to explore how geometry impacts hemodynamics, finding
that small changes in radius and length have significant impacts on
hemodynamics. As patient-specific modeling becomes increasingly popular in
medical applications, we argue that quantifying uncertainty in this context is
critical for interpretation of predictions.

Figure 1: Workflow for our methodology, consisting of image segmentation and
analysis, creation of a labeled directed tree, junction node modification,
change point analysis, radius identification, and construction of a
computational domain for fluid dynamics simulations.

##  2 Methods

###  2.1 Image analysis

We analyze two image types for this study, a chest CT image and an MRA of the
chest and neck region. The CT image is from a healthy, 67-year-old female
volunteer and captured using Omnipaque 350 contrast agent [37]. It can be
found in the publicly available Vascular Model Repository (see
http://simvascular.github.io/). The MRA is from a DORV patient from the Texas
Children’s Hospital Heart Center. Data collection for this study was approved
by the Baylor College of Medicine Institutional Review Board (H-46224: “Four-
Dimensional Flow Cardiovascular Magnetic Resonance for the Assessment of
Aortic Arch Properties in Single Ventricle Patients”). Time-resolved, contrast
enhanced MRA was performed using localizing sequences from the cardiac MRI
with a bolus injection of 0.1 mL/kg of intravenous gadolinium contrast. More
details can be found in [15].

3D rendering. We segment and reconstruct the pulmonary arterial and aortic
networks from their respective sources using the open-source segmentation
software 3D Slicer developed by Kitware, Inc. (see http://www.slicer.org) [38,
39]. We generate five segmentations of the pulmonary arteries and ten of the
aortic vessels, each from the same image using the same settings. The
segmentation is a 3D-rendered volume obtained using built-in tools from 3D
Slicer, including painting, thresholding, cutting, and islanding (Figure 2).
Image intensities in CT images range from $50-3027$ Hounsfeld units and
$75-327$ Hounsfeld units for the aortic MRA. From 3D Slicer, the segmentation
is saved to an STL format and opened in Paraview [40]. It is then converted to
a VTK polygonal data file to be processed through VMTK (see
http://www.vmtk.org/) [41].

Centerlines. We utilize VMTK [26] to determine centerlines in 3D-renderings of
the pulmonary arteries and aorta. This software uses maximally inscribed
spheres to determine centerlines in a network of vessels [41]. Spheres, which
touch at least four points on the 3D-rendering, are inscribed throughout each
vessel and the medial axis is approximated by a Voronoi diagram [41].
Centerlines are then defined as minimal paths along the inverse of the
spheres’ radii [41]. User-specified inlet and outlet points establish
centerline boundaries, guiding VMTK in a recursive process that works from
terminal to proximal vessels to define centerlines. When two centerlines
intersect, a junction point is placed. Centerlines within our segmentations
are shown in Figure 2. In many vessels, junction points defined by VMTK do not
occur in the center of vessel bifurcations, often occurring too early, leading
to overestimations in radii and underestimations in length.

###  2.2 Labeled directed tree

Large vessels. Existing algorithms were used in MATLAB to extract information
from VMTK and construct a labeled directed tree [11, 42]. Centerlines are
transformed to edges with $x,y,\text{ and }z$ coordinates and junction points
are denoted as nodes. A connectivity matrix, describing vessel radii, length,
and position, is constructed from the network topology. Radii are calculated
from maximally inscribed spheres in VMTK. Vessel length is calculated by
determining the Euclidean distances between $x,y,\text{ and }z$

Figure 2: 3D renderings from (a) normotensive pulmonary arterial network and
(b) the aortic vasculature of a DORV patient. Centerlines, shown as black
continuous lines, and structured trees constitute a labeled directed tree.

Table 1: Parameters used in the fluid dynamics model for the pulmonary and aortic networks. Quantity | Pulmonary | Aorta
---|---|---
$T$ (s) | $0.85$ | $0.8$
$\rho$ (g/cm3) | $1.055$ | $1.055$
$\mu_{L}$ (g/cm/s) | $0.032$ | $0.032$
$\mu_{S}(r_{0})$ | Equation (20) | Equation (20)
$k^{L}_{1}$ (g/cm/s2) | $2.5\times 10^{6}$ | $5.0\times 10^{6}$
$k^{L}_{2}$ (1/cm) | $-15$ | $-25$
$k^{L}_{3}$ (g/cm/s2) | $6.4\times 10^{4}$ | $6.0\times 10^{5}$
$k^{S}_{1}$ (g/cm/s2) | $2.5\times 10^{7}$ | $5.0\times 10^{6}$
$k^{S}_{2}$ (1/cm) | $-15$ | $-20$
$k^{S}_{3}$ (g/cm/s2) | $8.0\times 10^{5}$ | $1.0\times 10^{5}$
$r_{\text{min}}$ (cm) | $0.001$ | $0.001$
$(\alpha,\beta)$ | $(0.88,0.697)$ | $(0.90,0.60)$
$\ell_{rr}$ | $15.75$ | $50$

coordinates of each vessel. Terminal vessels are identified as degree-1
vessels with no further branching. The ten aortic networks each have $13$
manually identified vessels. In contrast, the rapidly branching structure and
wide range in spatial scale of the pulmonary arteries make it difficult to
extract a consistent number of vessels during each segmentation. The five
segmentations of the pulmonary arteries have between $167$ and $238$ vessels.
To enable direct comparison, we prune these networks to generate trees with
identical number of imaged vessels. We utilize a radius pruning method based
on [43], iteratively removing the smallest terminal vessels with terminal
siblings until each network reaches the desired number of vessels.

Small vessels. Due to limits in image resolution, it is difficult to obtain
measurements for small arteries and arterioles. To maintain the morphometric
accuracy of the system, we use asymmetrically-branching, self-similar
structured trees to extend the networks [44]. In the structured trees,
daughter vessel radii $r_{(d_{1},d_{2})}$ and lengths $l_{(d_{1},d_{2})}$ are
scaled relative to their parent $r_{p}$ such that

$\displaystyle r_{d_{1}}=\alpha r_{p},\ \ r_{d_{2}}=\beta r_{p},\ \
l_{(d_{1},d_{2})}=r_{(d_{1},d_{2})}\ell_{rr}$ (1)

where $\alpha,\beta<1,\alpha+\beta\geq 1$ and $\ell_{rr}$ is a length to
radius ratio. These values can be found in Table 1 and were based on
literature [44, 45]. The structured trees bifurcate until a specified minimum
radius, $r_{\text{min}}$, is reached. For this study, $r_{\text{min}}=0.001$
cm, which is consistent with a red blood cell’s diameter [46].

###  2.3 Junction node modification

We modify VMTK-derived nodes using an automated process that places them in a
physiologically accurate location. Our algorithm works recursively, starting
at terminal vessels and marching towards proximal vessels, systematically
reallocating centerline information. We first find the normalized distance
between each $x,y,$ and $z$ coordinate point along centerlines of the daughter
vessels. We isolate portions in daughter vessels that are within $10\%$ of
each other and place the new node at the end of this derived appendage.
Between the original and new node, we compute a pointwise average of the
daughter vessel coordinates and append this new line to the end of the parent
(Figure 3).

Next, we refine centerlines to ensure a gradual transition between parent and
daughter vessels, avoiding abrupt and angular adjustments. To smooth the
centerlines, we compute a weighted average in sections where daughters are
within $10-20\%$ of each other and reassign centerlines based on this average.
The degree of weighting is determined by the ratio of radius and length
between parent and daughter vessels. In the pulmonary arteries, we compute a
double-weighted average for each daughter in the system, accounting for
asymmetrical branching patterns. Since the aorta consists of a large, curved,
main branch, centerlines must lie closer to the main branch than the neck
vessels. We accomplish this by calculating a 100-times weighted average
towards the main branch. Based on these adjustments, new lengths and radii are
assigned to the labeled directed tree (Figure 3).

Figure 3: Junction node adjustment in 3D space. Black lines and teal circles
represent the centerlines and junctions placed by VMTK for (a) the pulmonary
vasculature and (b) the aorta. Pink lines and purple circles represent the
centerlines and junctions after our adjustment algorithm.

###  2.4 Statistical change point algorithm

In the context of determining the radius in each blood vessel, change points
can be used to determine locations along the vessel where there are notable
alterations in the radius. These alterations are induced by the non-smooth
segmentations, rapidly branching networks, and imprecise junction locations.
We utilize change points to identify vessel regions that best represent its
radius, avoiding extreme and uncertain sections.

In this study, a change point, $\psi\in\\{1,\dots,n-1\\}$, is identified
within an ordered sequence of radii $\mathbf{r}=\\{r_{1},\dots,r_{n}\\}$ along
the vessel for which statistical properties of $\\{r_{1},\dots,r_{\psi}\\}$
and $\\{r_{\psi+1},\dots,r_{n}\\}$ are significantly different. A dataset may
contain multiple change points, represented as
$\boldsymbol{\psi}=\\{\psi_{1},\psi_{2},\dots,\psi_{m}\\}$, where $m$ is the
number of change points. We employ a two-step process using publicly available
R packages to define the (i)number and (ii) location of change points in each
vessel. This approach allows us to create an automated framework, effective
for analyzing and comparing large datasets.

Number of change points. We use the R/changepoint package to determine the
optimal number of change points in each vessel [47]. This program employs
binary segmentation, iteratively partitioning the dataset based on fitting
piecewise constant models in each segment.

For each iteration, hypothesis testing is used to determine if a change point
should be placed within a segment. The null hypothesis, $H_{0}\\!:w=0$,
corresponds to no change points, and the alternative, $H_{a}\\!:w=1$, to one
change point [47]. To determine if a change has occurred, the maximum log-
likelihood is calculated as

$\displaystyle\lambda=2\left(\max_{\psi_{1}}M\\!L(\psi_{1})-\log\left(p(\mathbf{r}_{1:n}|\hat{\theta})\right)\right),$
(2)

where $\hat{\theta}$ is the maximum likelihood estimate of the parameters,
including the slope, variance, and mean [47], and the maximum likelihood
$M\\!L$ is

$\displaystyle
M\\!L(\psi_{1})=\log\\!\left(\\!p(\mathbf{r}_{1:\psi_{1}}|\hat{\theta}_{1})\right)+\log\\!\left(p(\mathbf{r}_{(\psi_{1}+1):n}|\hat{\theta}_{2})\right),$
(3)

where $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are parameter estimates before
and after the change point. The null hypothesis is rejected when $\lambda$ is
greater than a specified threshold. Upon rejection of the null hypothesis, a
change point is placed at $\max_{\psi_{1}}\\!M\\!L(\psi_{1})$ [47]. In this
study, we limit the number of change points to three to avoid overfitting and
enable consistent analysis between datasets.

Binary segmentation in R/changepoints efficiently identifies the number of
change points but is sensitive to initial partitioning, offering a low
computational cost at the expense of providing precise change point locations
[48, 49, 47]. Therefore, we employ a second software to determine change point
location, described below.

Change point placement. We use the customizable and flexible R/segmented
package to find optimal change point locations, maximizing the likelihood
between a linear regression model and the data[50, 51] (Figure 4). This
package requires that the user inputs the number of change points, obtained as
described above, and a regression model.

This study uses a linear regression model within R/segmented, which fits
piecewise linear functions between iteratively placed change points where
statistical properties, slope, mean, and variance, differ at a rate exceeding
a predetermined threshold. Iteration continues until there are no changes in
the maximum likelihood comparing the radii data and piecewise linear function.

The regression model between mean response, $E[\mathbf{r}]$, and variable
$\theta$ is defined as

$\displaystyle g(E[\mathbf{r}])=\eta(l)+b\times h(\hat{\theta};\psi_{i}),$ (4)

where $g(\cdot)$ is the link function for $E[\mathbf{r}]$,
$h(\hat{\theta};\psi_{i})$ is the non-linear parameterization for $\theta$,
$\psi_{i}$ is a change point with $i=1,2,3$, and $\eta(l)$ is the predictor. A
first-order Taylor expansion around an initial $\psi_{i}^{(0)}$ is used to
approximate

$\displaystyle h(\hat{\theta};\psi_{i})\approx
h(\hat{\theta};\psi_{i}^{(0)})+\left(\psi_{i}-\psi_{i}^{(0)}\right)h^{\prime}(\hat{\theta};\psi_{i}^{(0)}).$
(5)

The response and variable form a piecewise relationship with the terms

$\displaystyle a\theta+b(\theta-\psi_{i})_{+}.$ (6)

Note that $(\theta-\psi_{i})_{+}=(\theta-\psi_{i})\times(I(\theta>\psi_{i}))$
where $I(\cdot)$ is the indicator function. $I(\cdot)=1$ when the statement is
true, and $I(\cdot)=0$, otherwise. Using this parameterization,

  * •

$a$ is the slope of the line segment prior to the change point.

  * •

$b$ is the difference in slopes of the lines before and after the change
point.

  * •

$(a+b)$ is the slope of the line segment following the change point.

The change point is fixed at $U^{(s)}=(\theta-\psi_{i}^{(s)})_{+}$ and
$V^{(s)}=-I(\theta)>\psi_{i}^{(s)}$ and iterated $s$ times until there is no
improvement in the $M\\!L$ estimate. Using radii, $\mathbf{r}$ (cm), as the
response, and length, $l$ (cm), as the predictor, the linear regression model
is given by

$\displaystyle\mathbf{r}(l)=\beta_{0}-\beta_{1}l$ (7)

where $\beta_{0},\beta_{1}$ are linear coefficients. This form was chosen to
account for vessel tapering [52].

Figure 4: Demonstration of the statistical change point calculation placing up
to three change points.

###  2.5 Radius extraction

Using change point positions and piecewise linear fits between them, we detect
an optimal region, $\mathbf{r}_{\text{opt}}$, within each vessel’s data that
best represents its radius (Algorithm 1). In non-tapering vessels (i.e., the
pulmonary vasculature), we apply a single representative radius $\hat{r}$ to
the whole vessel by averaging over $\mathbf{r}_{\text{opt}}$. In tapering
vessels (i.e., the aortic vasculature), we assign the inlet,
$\hat{r}_{\text{in}}$, and outlet, $\hat{r}_{\text{out}}$, radii from an
exponentially decaying line fit to $\mathbf{r}_{\text{opt}}$. In both vessel
types, we designate the segment prior to the first change point as the
junction region, excluding it from the radius extraction process to avoid
outliers that could potentially distort radius choice. In the following
sections and Algorithm 1, we describe how $\mathbf{r}_{\text{opt}}$ is
selected. Figure 5 and the supplemental material provide additional pictorial
examples of this process.

Non-tapering vessels. To determine $\mathbf{r}_{\text{opt}}$, we compare the
slopes of regression lines between change points, choosing the section with
the smallest slope, $slp_{\min}$. For most vessels, this section is a stable
representation of a vessel’s features as it is minimally affected by
oscillations in data caused by non-smooth segmentations, overestimation within
the junction region, and underestimation from image resolution limits. If the
data encompassed by $slp_{\min}$ accounts for more than $25\%$ of observed
vessel radius measurements, we consider this to representative of the vessel’s
radius, assigning it to be $\mathbf{r}_{\text{opt}}$. Enforcing this condition
excludes small flat regions that may be present near the beginning or end of
the vessel.

If $slp_{\min}$ includes less than $25\%$ of data points, we adopt an
alternative approach to find $\mathbf{r}_{\text{opt}}$. For vessels where the
first change point, $\psi_{1}$, is positioned beyond the $25\%$ mark of the
vessel, $\mathbf{r}_{\text{opt}}$ begins here. If $\psi_{1}$ is located before
the $25\%$ mark, the first change point may be in the junction region, so we
start $\mathbf{r}_{\text{opt}}$ at the second change point, $\psi_{2}$. To
avoid excluding data, when $\psi_{2}$ is located more than $50\%$ into the
vessel, we instead start $\mathbf{r}_{\text{opt}}$ at the $25\%$ mark in the
vessel.

Next, our algorithm chooses the end of $\mathbf{r}_{\text{opt}}$, considering
either the last change point, $\psi_{m}$, or the last radius observation,
$r_{n}$, depending on the data structure. Due to imperfections in the rendered
image, such as noise and poor contrast, the ends of some vessels, particularly
terminal vessels, may not be accurately captured [53, 54]. As a result, radius
observations after $\psi_{m}$ often differ significantly from the rest of the
data. To account for this, we compare the slopes in the second-to-last and
last sections of data, $slp_{m}$ and $slp_{m+1}$. We identify vessels for
which $slp_{m}$ and $slp_{m+1}$ are significantly different by determining the
relative error between the two as

$R\\!E=\left|\frac{slp_{m+1}-slp_{m}}{slp_{m}}\right|.$ (8)

If $R\\!E$ is greater than or equal to a threshold, determined by the data, we
omit radius observations after $\psi_{m}$, ending $\mathbf{r}_{\text{opt}}$ at
$\psi_{m}$. Otherwise, $r_{n}$ is chosen as the end point.

There are several vessels that do not follow the aforementioned patterns and
require separate consideration. For vessels with no detected change points, we
compute the radius by averaging over the entire vessel. The root vessel of the
pulmonary arteries also has distinct characteristics. Manual selection of
inlets in VMTK results in qualitatively different data for this vessel.
Therefore, the radius of the main pulmonary artery (MPA) is determined from
the center $50\%$ of the vessel.

Several studies of vascular morphometry have found that for small vessels in
healthy subjects, the radius of a parent vessel is greater than each of its
daughters and the area of the daughters is greater than their parent [55, 56,
57, 46]

$\displaystyle r_{p}$ $\displaystyle>r_{d_{1}},r_{d_{2}},\ \ \ \
A_{p}<A_{d_{1}}+A_{d_{2}}.$ (9)

In vessels where these conditions are not satisfied, we adjust radii depending
on the vessel’s position within the network. When a terminal vessel violates
Equation (9), we assign a new radius value using the structured tree framework
outlined in Section 2.2. If a single terminal daughter, $d_{1}$, violates the
assumption, we assign it new radius of $r_{d_{1}}=\alpha r_{p}$. If both
daughters violate the assumption, we assign new radii values of
$r_{d_{1}}=\alpha r_{p}$ and $r_{d_{2}}=\beta r_{p}$. Internal vessels that
violate Equation (9) are manually assigned a new radius by choosing a more
suitable value from inspection of the data.

Tapering vessels. The methodology for determining $\mathbf{r}_{\text{opt}}$
for tapering vessels is similar. Instead of using the minimum slope, to
determine the optimal region, we use the longest interval between change
points, $l_{\max}=\max(\mathbf{sec})$, where
$\mathbf{\mathbf{sec}}=\\{l_{1},\dots,l_{M+1}\\}$ represents sections of data
between change points and $l_{1:M+1}$ represents the length of each section.
This region must include over $40\%$ of radii observations in a given vessel;
otherwise, we follow the same process as non-tapering vessels.

The equation for exponentially tapering vessels is

$\mathbf{r}(l)=n_{1}\exp(-n_{2}l)+n_{3},$ (10)

where $n_{1}$ is $r_{\text{in}}-r_{\text{out}}$ $r_{\text{in}}$, $n_{2}$ is
the average degree of taper in a vessel, and $n_{3}$ is the outlet radius
[15]. $n_{1,2,3}$ are optimized using MATLAB’s function fminsearch() [58] to
fit an exponentially tapering function along the length of the vessel,
minimizing the maximum squared error over$\mathbf{r}_{\text{opt}}$. Values for
$\hat{r}_{\text{in}}$ and $\hat{r}_{\text{out}}$ are determined by predicted
radius values from the exponential fit at the first point and last point of
the optimal region, respectively.

###  2.6 Fluid dynamics model

We predict hemodynamics in each vessel using a 1D fluid dynamics model that
relates blood pressure, blood flow, and cross-sectional area. The model is
separated into two domains predicting hemodynamics in: (1) large vessels by
solving the cylindrical Navier$-$Stokes equations in a reconstructed geometry
from imaging data, and (2) small vessels by solving a linearized model in
structured trees. Parameters for each domain are listed in Table 1 and the
supplement.

Large vessels. Similar to [15, 16, 44, 11], we predict blood pressure $p(x,t)$
(mmHg), blood flow $q(x,t)$ (cm3/s), and area deformation $A(x,t)$ (cm2) in
large vessels by solving the 1D Navier$-$Stokes equations. We assume that each
vessel can be represented by a cylindrical and axisymmetric tube with a
circular cross-section and impermeable wall. We assume that blood is
incompressible, Newtonian, viscous, and homogeneous and that the flow is
irrotational and laminar.

For this system, mass conservation and momentum balance are given by

$\displaystyle\frac{\partial A}{\partial t}+\frac{\partial q}{\partial x}$
$\displaystyle=0$ (11) $\displaystyle\frac{\partial q}{\partial
t}+\frac{\partial}{\partial
x}\left(\frac{q^{2}}{A}\right)+\frac{A}{\rho}\frac{\partial p}{\partial x}$
$\displaystyle=\frac{-2\pi\nu R}{\delta}\frac{q}{A},$ (12)

where $\mu_{L}$ is viscosity, $\nu=\mu_{L}/\rho$ is kinematic viscosity, $R$
(cm) is the radius, $t$ (s) denotes the temporal coordinate, and $x$ (cm) is
the axial position. We assume a flat velocity profile with a linearly
decreasing boundary layer [59, 60, 61, 44] with thickness $\delta=\sqrt{\nu
T/2\pi}$ (cm) where $T$ (s) is the length of the cardiac cycle, such that

$\displaystyle u_{x}(r,x,t)=\begin{cases}\bar{u}_{x},\ \ \ r<R-\delta\\\
\frac{\bar{u}_{x}}{\delta}(R-r),\ \ \ R-\delta<r\leq R.\end{cases}$ (13)

1 Function  _get_region_ :

input : $\textbf{r}=\\{r_{1},\dots,r_{n}\\}$,
$\boldsymbol{\psi}=\\{\psi_{1},\dots,\psi_{m}\\},\textbf{slp}=\\{slp_{1},\dots,slp_{m+1}\\},\textbf{sec}=\\{l_{1},\dots,l_{m+1}\\}$

output : $\textbf{r}_{\text{opt}}$

2

3 if  _vessels do not taper_ then

4 if  _$m=0$_ then

5 $\textbf{r}_{\text{opt}}=\textbf{r}$

6 else

7 $slp_{\min}=\min\left(\\{|slp_{2}|,\dots,|slp_{m+1}|\\}\right)$

8 $n_{s}$ = length($slp_{\min})$

9 $k=\text{index}(\psi_{\min-1})$

10 if  _$n_{s}/n\geq 0.25$_ then

11 $\textbf{r}_{\text{opt}}=\textbf{r}_{k:k+n_{s}}$

12 else if  _index $(\psi_{1})/n\geq 0.25$ _ then

13
$\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\text{index}(\psi_{1}),\text{index}(\psi_{m}))$

14 else if  _index $(\psi_{2})/n\geq 0.5$_ then

15 $\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\lceil
n/4\rceil,\text{index}(\psi_{m}))$

16 else

17
$\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\text{index}(\psi_{2}),\text{index}(\psi_{m}))$

18 end if

19

20 end if

21

22 else if  _vessels taper_ then

23 $l_{max}=\max\left(\text{sec}\right)$

24 $k$ = index$(\psi_{l_{max}-1})$

25 if  _$l_{max}/n\geq 0.4$_ then

26 $\textbf{r}_{\text{opt}}=\textbf{r}_{k:k+l_{max}}$

27 else if  _index $(\psi_{1})/n\geq 0.4$ _ then

28
$\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\text{index}(\psi_{1}),\text{index}(\psi_{m}))$

29 else if  _index $(\psi_{2})/n\geq 0.5$_ then

30 $\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\lceil
n/4\rceil,\text{index}(\psi_{m}))$

31 else

32
$\textbf{r}_{\text{opt}}=\text{slopes}(\textbf{r},\textbf{slp},\text{index}(\psi_{2}),\text{index}(\psi_{m}))$

33 end if

34

35 end if

36

37

38 Function  _slopes_ :

input : $\textbf{r}=\\{r_{1},r_{2},\dots,r_{n}\\}$,
$\textbf{slp}=\\{slp_{1},slp_{2},\dots,slp_{m+1}\\},i_{1},i_{2}$

output : $\textbf{r}_{\text{opt}}$

39

40 if  _$|(slp_{m+1}-slp_{m})/slp_{m}|\geq\text{threshold}$ _ then

41 $\textbf{r}_{\text{opt}}=\textbf{r}_{i_{1}:i_{2}}$

42 else

43 $\textbf{r}_{\text{opt}}=\textbf{r}_{i_{1}:n}$

44 end if

45

46

Algorithm 1 Given radius data $\mathbf{r}$, change points $\boldsymbol{\psi}$,
slopes of linear fits between change points $\mathbf{slp}$, and the length of
the sections between change points, $\mathbf{sec}$, these functions extract an
optimal region from radius data to determine a representative vessel radius,
$\mathbf{r}_{\text{opt}}$. In non-tapering vessels, $\mathbf{r}_{\text{opt}}$
is averaged and applied to the entire vessel. In tapering vessels, data from
the region is fit to an exponential curve.

To close the system of equations, we model the vessel wall as a linearly
elastic and isotropic material, allowing us to relate pressure and area [62]

$\displaystyle
p(x,t)-p_{0}=\frac{4}{3}\frac{Eh}{r_{0}}\left(\sqrt{\frac{A_{0}}{A}}-1\right)$
(14)

where $E$ (g/cm/s2) is Young’s modulus, $h$ (cm) is vessel wall thickness,
$p_{0}$ (mmHg) is a reference pressure, $A_{0}$ (cm2) is the corresponding
cross-sectional area, and $R_{0}$ (cm) is the radius of the vessel determined
by our extraction algorithm. To account for changes in wall composition,
vascular stiffness is modeled as exponentially increasing with vessel radius
such that

$\displaystyle\frac{Eh}{r_{0}}=k^{L}_{1}e^{k^{L}_{2}r_{0}}+k^{L}_{3},$ (15)

where $k^{L}_{1}$ (g/cm/s2), $k^{L}_{2}$ (1/cm), and $k^{L}_{3}$ (g/cm/s2) are
constants based on literature [42, 63, 15].

Figure 5: Radius extraction for (a) pulmonary arteries and (b) aortic vessels.
Change points (red dots) are placed within a vessel’s radius data and we
determine a representative radius for each vessel (green) by analyzing the
data’s structure, avoiding the junction region (pink box). In the aorta,
estimated $r_{\text{in}}$ and $r_{\text{out}}$ are denoted by blue stars. We
find the standard deviation of data in the optimal region and construct a CDF
and PDF, justifying a normal distribution fit.

Since the system of equations is hyperbolic, we define boundary conditions at
each end of the vessel. At the inlet of the root vessel, we prescribe an
inflow waveform from data. In the pulmonary arteries, we specify a flow
waveform at the MPA using control MRI data available from SimVascular [28]. In
the aortic vasculature, a flow waveform at the root of the ascending aorta is
extracted from 4D-MRI data provided by Baylor College of Medicine and Texas
Children’s Hospital. At junctions, we enforce continuity of pressure and
conservation of flow between parents, $p$, and daughters, $d_{1}$ and $d_{2}$

$\displaystyle p_{p}(L,t)$ $\displaystyle=p_{d_{1}}(0,t)=p_{d_{2}}(0,t)$ (16)
$\displaystyle q_{p}(L,t)$ $\displaystyle=q_{d_{1}}(0,t)+q_{d_{2}}(0,t).$ (17)

At the end of each terminal large vessel, we use a structured tree to match
impedance values. This model is solved using the two-step Lax-Wendroff finite
difference scheme [11, 42, 16, 15, 44, 45].

Small vessels. In small vessels, viscous forces are dominant, allowing us to
neglect the nonlinear inertial terms and linearize Equations (11) and (12)
such that

$\displaystyle i\omega Q+\frac{A_{0}(1-F_{J})}{\rho}\frac{\partial P}{\partial
x}$ $\displaystyle=0,F_{J}=\frac{2J_{1}(w_{0})}{w_{0}J_{0}(w_{0})}$ (18)
$\displaystyle i\omega CP+\frac{\partial Q}{\partial x}$
$\displaystyle=0,C=\frac{\partial p}{\partial
A}\approx\frac{3}{2}\frac{r_{0}}{Eh}$ (19)

where $J_{0},J_{1}$ are the zeroth and first order Bessel functions,
$w_{0}^{2}=i^{3}r_{0}\omega/\nu$, $C$ is vessel compliance, and $Eh/r_{0}$ has
the same form as Equation (15) but with parameters $k^{S}_{1},k^{S}_{2},$ and
$k^{S}_{3}$. Viscosity becomes more significant for smaller vessels [61], so
in this domain we compute viscosity as a function of vessel radius

$\displaystyle\mu_{S}(r_{0})$
$\displaystyle=\frac{D\mu_{L}}{3.2}\left(1+D(\mu_{0.45}-1)\right)$ (20)
$\displaystyle\mu_{0.45}(r_{0})$
$\displaystyle=6e^{-0.17r_{0}}+3.2-2.44e^{-0.12r_{0}^{0.645}},$ (21)

where $D=\left(2r_{0}/(2r_{0}-1.1)\right)^{2}$, and $\mu_{0.45}(r_{0})$ is the
relative viscosity at a hematocrit level of 0.45.

The linearized equations are converted to discrete, frequency-domain Fourier
series, which results in a wave equation

$\displaystyle\frac{\omega^{2}}{c}Q+\frac{\partial^{2}Q}{\partial
x^{2}}=0,c=\sqrt{\frac{A_{0}(1-F_{J})}{\rho C}}$ (22)

that can be solved analytically for flow $Q$ and pressure $P$

$\displaystyle Q(x,\omega)$ $\displaystyle=a\cos(\omega x/c)+b\sin(\omega
x/c)$ (23) $\displaystyle P(x,\omega)$
$\displaystyle=\frac{i}{g_{\omega}}\left(-a\sin(\omega x/c)+b\cos(\omega
x/c)\right),$ (24)

where $a$ and $b$ are integration constants and
$g_{\omega}=\sqrt{CA_{0}(1-F_{J})/\rho}$ [15, 42, 44]. We relate $Q$ and $P$
through impedance, $Z(x,\omega)=P(x,\omega)/Q(x,\omega)$, with the input
impedance at $x=0$

$\displaystyle Z(0,\omega)=\frac{ig_{\omega}^{-1}\sin(\omega
L/c)+Z(L,\omega)\cos(\omega L/c)}{\cos(\omega
L/c)+ig_{\omega}Z(L,\omega)\sin(\omega L/c)}$ (25)

for $\omega\neq 0$ and $x=L$ representing the end of the vessel. The zero-
frequency value is defined as

$\displaystyle\lim_{\omega\rightarrow
0}Z(0,\omega)=\frac{8\mu_{S}\ell_{rr}}{\pi r_{0}^{3}}+Z(L,0).$ (26)

For boundary conditions, we enforce pressure continuity and flow conservation
at bifurcations such that

$\displaystyle\frac{1}{Z_{p}}=\frac{1}{Z_{d_{1}}}+\frac{1}{Z_{d_{2}}},$ (27)

where $Z_{p},Z_{d_{1}}$, and $Z_{d_{2}}$ are the impedance values of the
parent and daughter vessels, respectively. We assume the terminal impedance at
the end of the structured tree is zero and use this to back-solve the
impedance values for all vessels in the system.

###  2.7 Uncertainty quantification

Statistical analysis. To understand image segmentation uncertainty, we perform
a one-way analysis of variance (ANOVA) test on five segmentations of the
pulmonary vasculature and ten segmentations of aortic vasculature before
modifying junctions or extracting radii based on change points. This
statistical test determines if there are significant differences among the
variance of different samples [64]. We conduct ANOVA tests in a vessel-
specific manner to determine if there are statistically significant
differences in vessel radii and length from different 3D-renderings of the
same subject. One-way ANOVA calculations are performed in R using the function
anova() from the R/stats package [65] (Table 2). More details about ANOVA can
be found in the supplemental material (S4).

To examine variability within the data, we calculate coefficients of variation
as CV$=\sigma_{\hat{r}}/\hat{r}$ for vessel radii in each segmentation [66]
(Figure 6). Here, $\hat{r}$ is the mean of the optimal radius region extracted
in Section 2.5 and $\sigma_{\hat{r}}$ is its standard deviation.

We also compare kernel density estimations (KDE) of vessel radii for each
segmentation to view network variability. The kernel density represents the
data using the probability density function (PDF) of the radii. These KDEs are
plotted to visualize the estimated distribution of the radii throughout the
segmented vasculatures [67] (Figure 6). We use MATLAB’s ksdensity() to convert
the radii from each segmentation into a PDF [58].

Fluid dynamics. We investigate uncertainty in fluid dynamics predictions by
varying image segmentation, junction placement, and vessel radius. To explore
the impact of segmentation, we predict hemodynamics in five distinct
segmentations of the pulmonary arteries and ten of the aorta from the same
subject. The complexity and small spatial scale of the pulmonary arteries
leads to variability in the vessel count in each segmentation. To enable
direct comparison between segmentations, we apply a radius pruning procedure
as described in Section 2.2, elucidating the impact of network size (Figure
7). We also apply our junction modification method and compare hemodynamics
within and between segmentations (Figure 8).

We select a representative segmentation to study the influence of vessel
radius on hemodynamics. To explore how radius selection propagates to fluid
dynamics, we simulate various radii choices in our representative networks. We
fit a normal distribution $\mathcal{N}(\hat{r},\sigma_{\hat{r}}^{2})$ to the
region selected in Section 2.5 and sample independently for each vessel
(Figure 5). Sampling from a normal distribution is justified by systematically
viewing PDF plots of radius data in the optimal region extracted. In vessels
with reassigned radii due to violating Equation (9), we sample from the
distribution $\mathcal{N}(\kappa\hat{r}_{p},\kappa\sigma_{\hat{r}_{p}}^{2})$,
where $\kappa=\frac{\hat{r}_{d}}{\hat{r}_{p}}$ is a parent-daughter scaling
factor. We sample $1000$ times in representative networks (Figure 9) and $100$
times in other segmentations (see supplemental material).

##  3 Results

Statistical analysis. One way ANOVA results (Table 2) reveal that there are
statistically significant differences in vessels between segmentations. This
means that the variability observed between each segmentation is larger than
what would be expected by solely random chance [64]. The ANOVA analysis for
seven pulmonary arteries matched in each segmentation show statistically
significant differences in length and radii. In the aortic network, the radii
in the innominate, brachial, right common carotid (RCC), and right subclavian
differ significantly and the length of the left common carotid (LCC),
brachial, left vertebral, and RCC also differ significantly between
segmentations. This analysis demonstrates that the raw data obtained from
3D-renderings is uncertain.

Table 2: One way ANOVA test results for error sum of squares (SSE) and
p-values for seven vessels in the pulmonary vasculature and eight vessels from
the aortic network that can be matched in each segmentation. For continuity,
all aortic vessels (ascending, arch I, arch II, and descending) were treated
as one exponentially tapering vessel. P-values in bold text are statistically
significant under a 0.05 threshold.

| Pulmonary | Aortic
---|---|---
| Name | SSE | p-value | Name | SSE | p-value
Radius ANOVA results | MPA | $2.82$ | $\mathbf{<0.001}$ | Aorta | $0.94$ | $0.11$
RPA | $5.45$ | $\mathbf{<0.001}$ | Inomminate | $2.73$ | $\mathbf{<0.001}$
LPA | $3.55$ | $\mathbf{<0.001}$ | LCC | $0.14$ | $0.97$
RIA | $3.20$ | $\mathbf{<0.001}$ | L Subclavian | $0.16$ | $0.32$
LIA | $13.75$ | $\mathbf{<0.001}$ | L Brachial | $0.10$ | $\mathbf{0.0093}$
RTA | $1.05$ | $\mathbf{<0.001}$ | L Vertebral | $0.004$ | $0.10$
LTA | $18.63$ | $\mathbf{<0.001}$ | RCC | $0.12$ | $\mathbf{<0.001}$
|  |  | R Subclavian | $0.12$ | $\mathbf{<0.001}$
Length ANOVA results | MPA | $246.59$ | $\mathbf{<0.001}$ | Aorta | $652$ | $0.40$
RPA | $177.56$ | $\mathbf{<0.001}$ | Inomminate | $0.77$ | $0.99$
LPA | $73.96$ | $\mathbf{<0.001}$ | LCC | $298$ | $\mathbf{0.02}$
RIA | $37.83$ | $\mathbf{<0.001}$ | L Subclavian | $4$ | $0.99$
LIA | $195.31$ | $\mathbf{<0.001}$ | L Brachial | $43.10$ | $\mathbf{0.047}$
RTA | $14.28$ | $\mathbf{<0.001}$ | L Vertebral | $58.59$ | $\mathbf{<0.001}$
LTA | $32.01$ | $\mathbf{<0.001}$ | RCC | $371$ | $\mathbf{<0.001}$
|  |  | R Subclavian | $9.37$ | $0.83$

Main pulmonary artery (MPA), right (RPA) and left (LPA) pulmonary arteries,
right (RIA) and left (LIA) interlobular arteries, and right (RTA) and left
(LTA) trunk arteries.

Figure 6 shows the CV and KDE for the pulmonary and aortic vasculatures. The
largest vessels in the pulmonary vasculature, with radii between $1$ and $1.5$
cm, have relatively small CVs, indicating a more certain measurement from
imaging data. For vessels with radii less than $0.5$ cm, CVs range from
$0-0.6$ with one outlier reaching a CV of over $1$. This vessel has a standard
deviation greater than its mean, which can occur in extreme cases where
terminal vessels are close to image resolution limits. In each segmentation,
smaller vessels consistently have higher CV values and a broader range of CVs.
In the aortic vasculature, CVs range from $0-0.32$. Vessels with radii less
than $0.4$ cm, and greater than $1$ cm exhibit CV values of less than $0.18$,
but vessels with radii between $0.4-1.0$ cm have CVs greater than $0.27$.
Similar to the pulmonary arteries, the largest aortic vessels have the
smallest CV of less than $0.1$.

The pulmonary KDE plots show a unimodal, right-skewed distribution with the
majority of vessel radii from all segmentations being less than $0.5$ cm. In
each segmentation, most vessels have radii less than $0.5$ cm peaking around
$0.1$ cm. However, the exact distribution of radii values varies across
networks, even after pruning. For example, in one segmentation, a radius of
$0.1$ cm is observed in $17.5\%$ of vessels, while in another, it is found in
only $12.5\%$. The latter network demonstrates a slightly more distributed
spread of radii values throughout its vessels. For the aortic networks, the
KDE follows a right-skewed bimodal distribution with a major mode at $0.25$ cm
and minor mode at $0.6$ cm. Based on the peaks in KDE, most aortic vessels
have a radius of less than $0.5$ cm with another large portion having a radius
between $0.5-0.75$ cm. Each segmentation follows a right-skewed bimodal
pattern, but similar to the pulmonary vessels, distinct networks have
different radii distributions. In particular, two segmentations have peaks at
approximately $0.25$ and $0.6$ cm. For one segmentation, the frequency of
these modes are $10\%$ and $4.5\%$, and for another segmentation, the
frequencies are $5\%$ and $4\%$, respectively.

Network pruning. Figure 7 shows the impact of standardizing network size
between segmentations in pressure and flow at the midpoint of the main (MPA),
right (RPA), and left (LPA) pulmonary arteries. Prior to pruning, there is a
wide range in predicted pressures between each segmentation with systolic
differences of up to $15$ mmHg. The segmentation with the highest number of
large vessels pre-pruning has a pressure range of between $18-31$ mmHg. The
segmentation with the smallest number of large vessels has a pressure range of
$5-15$ mmHg. Once all networks were pruned, the pressure range between each
segmentation lessened to about $5$ mmHg. Networks with the most vessels pre-
pruning had more significant drops in pressure, likely because the total
resistance of the network follows the vessel count. There are still
differences between pressure values in distinct segmentations, but these
differences are less pronounced and can likely be attributed to segmentation
uncertainty in the smallest vessels. While pruning has significant impacts on
pressure distribution, there are less pronounced differences in volumetric
flow. Nevertheless, both segmentation and pruning impact flow dynamics. In the
RPA, before pruning there is a peak flow range of nearly $100$ mL/s and after
pruning, this range becomes less than $50$ mL/s. Discrepancies in flow in the
MPA and LPA are smaller.

Junction node modification. Pressure and flow waveforms calculated at vessel
midpoints before and after modifying junctions are shown in Figure 8. The
variation before and after modifying junctions is significantly smaller in the
pulmonary vessels compared to the aorta. In most cases there is a difference
of $1-3$ mmHg between pressure predictions in the MPA, RPA, and LPA. The
modification impacts pulse pressure more than systolic or diastolic pressure.
The difference between flow predictions is also smaller in the pulmonary than
the aortic vasculature, but there is a distinct effect on the shape of the
flow predictions. Before modifying junctions, pulse pressure in the aortic
vasculature is about $80$ mmHg with values ranging from $40-120$ mmHg. After
this change, the pulse pressure is $50$ mmHg, ranging from $50-110$ mmHg.
Before modifying junctions, there is increased flow to the branching vessels
(LCC and brachial artery), whereas after this modification, the flow increases
in the main aortic vessels (descending aorta). Differences in the reflection
waves between predictions before and junction modification are also observed.
In the LCC, before junction modification, there are reflections between
$0.4-0.5$ s and again at the end of the cardiac cycle. The solid black line
shows a large negative dip and then one large reflection around $0.6$ s. After
modifying junctions, we observe one steady reflection taking place between
$0.5-0.7$ s. The opposite trend is apparent in the brachial artery, where
predictions before junction modification show one reflection and those after
show two smaller reflections. Note that flow waveforms after modifying
junctions are a close to fit to this patient’s waveforms extracted from their
4D-MRI (see supplemental material, S5).

Change point analysis. Figures 4 and 5 depict change point analysis used to
identify vessel radii. We choose a maximum of three change points in each
vessel to balance model complexity and interpretability. We analyze over 1000
vessels within five segmentations of pulmonary arteries and ten of the aorta,
finding that a majority of vessels have an optimal number of three change
points. Approximately $9\%$ of pulmonary arteries and $33\%$ of aortic vessels
have less than three change points. Our change point and radius extraction
methods enable robust analysis of large datasets with less than $1\%$ of
vessels requiring post-processing to satisfy Equation 9.

Varying radii. Figure 9 shows how varying vessel radii impacts pressure and
flow waveforms in a representative network. Pulmonary pressure predictions in
the MPA, RPA, and LPA for sampled radii range from $10-25$ mmHg in systole. In
the MPA, there is slight variation in flow predictions in systole ranging from
approximately $275-290$ mL/s. In the RPA and LPA, variations are more
pronounced ranging from $100-175$ mL/s and $50-100$ mL/s, respectively. Aortic
pressures range from $90-160$ mmHg in systole. Flow predictions in the
descending aorta range from $250-350$ mL/s in stystole, while both the LCC and
brachial range from $30-40$ mL/s.

Figure 6: CV and KDE for (a) pulmonary arteries and (b) aortic vessels. Each
color represents a distinct segmentation. Figure 7: Pressure and flow dynamics
over a cardiac cycle at the midpoint of the MPA, RPA, and LPA in the pulmonary
arteries before and after standardizing network size through pruning. Each
color represents a distinct segmentation of the same network, solid lines show
hemodynamics prior to pruning, and dotted lines display results after pruning.
Figure 8: Pressure and flow dynamics over a cardiac cycle at the midpoint of
(a) the MPA, RPA, and LPA in the pulmonary arteries and (b) the descending
aorta, LCC, and brachial artery before and after modifying VMTK-derived
junctions. Each color represents a distinct segmentation of the same network,
solid lines show hemodynamics prior to modifying junctions, and dotted lines
display results after changing junctions. Figure 9: Pressure and flow dynamics
over a cardiac cycle at the midpoint of (a) the MPA, RPA, and LPA in the
pulmonary arteries and (b) the descending aorta, LCC, and brachial artery.
Gray lines (1000) denote simulations with radii sampled from a normal
distribution based on the mean and standard deviation of each vessel’s radius
and the colored line (1) denotes the simulations with radii set at the mean
radius based on our extraction procedure. These simulations were based on a
representative segmentation deemed of the pulmonary and aortic vasculatures.

##  4 Discussion

In this study, we devise novel algorithms for generating labeled directed
trees from medical images of the pulmonary arteries and aortic vasculature. We
develop a robust method for accurate placement of junctions in vessel
bifurcations and use change points to extract vessel radii for use in a CFD
model. To test our approach, we generate multiple segmentations of identical
vasculatures, implement a standardization procedure to generate comparable
networks, and explore the impact of different radii and length in each vessel.
We find that varying geometric features, 3D-renderings, and network size have
significant impacts on pressure and flow predictions in both the pulmonary and
aortic vessels.

Small vessel variability. Compared to aortic vessels, the pulmonary
vasculature has larger CVs, especially in the smallest vessels (Figure 6). A
high CV indicates that the vessel exhibits significant variability in its
radius along its length, signifying an irregular structure. This variability
can be due to physiological factors, such as genetic anomalies, disease,
inflammation, or local blood flow patterns [68, 69, 70, 71, 72, 73], and
measurement errors, due to limited image resolution, segmentation
inaccuracies, or inexactness in centerline generation [74, 7, 9, 10]. It is
well known that imaging variables, such as noise, poor contrast, and blurred
boundaries, increase the difficulty of extracting the smallest vessels [53,
54, 75]. In a study by Colebank et al. [11] looking at uncertainty in image
segmentation, CV increases as vessel radius decreases, suggesting that small
vessel dimensions are sensitive to segmentation parameters. Horssen et al.
[76] postulate that small vessel dimensions are most likely underestimated due
to limits in image resolution. Schwarz et al. [77] find that in the human
coronary arterial circulation, the smallest vessels are subject to
segmentation errors, resulting in triangular loops and spurious side branches
that need to be corrected in post-processing. These studies validate our
findings as vessels in the aortic network are much larger in both radius and
length than the rapidly branching pulmonary network. Despite imaging
constraints, our innovative approach for determining vessel dimensions
mitigates uncertainty because it systematically bypasses artifacts of the
imaging process.

The bimodal nature of the KDE plots in aortic vessels (Figure 6) is due to the
subset of aortic vessels (ascending aorta, arch I, arch II, and descending
aorta) having significantly larger radii than the neck vessels. In contrast,
the pulmonary vasculature has a rapidly branching structure, giving rise to a
unimodal right-skewed PDF. The importance of quantifying uncertainty in image
segmentation through the use of PDFs has been supported by several studies. In
Baumgartner et al. [78], they generate PDFs for each segmentation and compare
it to the distribution of the actual image. They argue that knowing these PDFs
can greatly improve image segmentation accuracy [78]. KDEs have also been used
for machine learning applications of image segmentations. Wang et al. [79]
leverage KDEs of abdominal blood vessels from CT angiograms to learn target
distributions and estimate vessel locations in segmentations of the same
vessels in different patients. The present study uses KDEs to understand
vessel distributions within and between segmentations and we note that this
analysis can be a powerful tool to gain information about vessel uncertainty.

Segmentation discrepancies. KDE and ANOVA results reveal large variation in
each segmentation of the pulmonary and aortic vasculatures (Figure 6, Table
2). We observe these discrepancies even though we control for 3D Slicer user,
segmentation subject, and threshold parameters. Our findings validate the
importance of investigating the segmentation and network generation process,
as they confirm complex geometries are susceptible to inconsistencies. KDEs
have been used to estimate the distribution of vessel radii in studies such as
[11, 80]. Even when using the same threshold, segmentations can vary due to
manual user error. In Mharib et al. [81] both semi-automatic and automatic
segmentation methods require user input, inviting potential errors. Besides
user error, most machine learning methods require vast amounts of training
data and performance decreases on large, complex networks due to resolution
constraints [81]. The results of our ANOVA tests also confirm there are
discrepancies between segmentations pre-junction modification. These show that
the radii and length of branching vessels have statistically significant
differences between segmentations. All pulmonary vessels compared exhibit
statistically significant differences between segmentations, whereas in the
aortic network, several vessels compared did not show statistically
significant differences. This is plausible because aortic vessels have larger
radii and length compared to pulmonary vessels, making it easier to generate
consistent segmentations.

Network standardization. Standardizing networks to have a consistent number of
vessels is imperative to enable comparison. After pruning pulmonary networks,
we find that differences between segmentations become less apparent. Once
vessels with an initially high branch count are pruned, their predicted
pressure decreases (Figure 7). This result likely stems from applying the same
boundary condition to each network assuming that continued branching follows a
structured network. One way to account for differences in the large vessel
network size is to adjust structured tree parameters, $\alpha$ and $\beta$,
generating networks with similar total vessel count. Miller et al. [43] argue
that networks with disparate branch numbers are an imaging consequence rather
than a biological factor and thus, pruning is essential when comparing
different segmentations. In addition, they find that more information about
networks can be gained by applying pruning techniques [43]. Colebank et al.
[11] find that choices in segmentation parameters greatly impact the number of
vessel in the network and that network size effects model output. Since mean
pressure and flow metrics are used to diagnose diseases, such as pulmonary
hypertension [82], standardizing results and assessing uncertainty is
critical. This also highlights the topic of determining the minimum number of
imaged arteries needed to explicitly model to obtain reliable diagnostic
results. If a subset of larger vessels can yield valid, comparable results, it
may be beneficial sense to exclude smaller vessels that exhibit drastic
imaged-related uncertainties. This idea has been studied by several groups,
who use surrogate or reduced order models to lessen computational complexity
[83, 84, 85, 86].

Junction node placement. The placement of nodes within junctions impacts both
length and radii of vessels. Large vessels, especially the aorta, are known to
taper along their length [52]. Due to the rapid change in vessel radii, node
placement in these vessels is increasingly important. Our algorithm addresses
this by shifting the VMTK-derived nodes down the vasculature, lengthening the
parent vessel and in turn, changing the radius data considered for each
vessel. Sections of centerline data that were previously characterized as part
of daughter vessels become reclassified as part of the parent. After this
reassignment, we pass the new radius data through our statistical change point
and radius extraction algorithms, which can change the radius used in
simulations depending on the data’s structure. We find that modifying nodes
impacts hemodynamics in both the pulmonary and aortic vasculature (Figure 8),
but it has a larger impact on aortic vessels. Since aortic vessels are longer
and wider than pulmonary vessels, geometric changes are more pronounced after
modifying junction placement and this effect is propagated to hemodynamics. In
particular, in the descending aorta, one segmentation has a variation of
nearly $60$ mmHg in systole between node locations. This is not the case in
pulmonary vessels, which have an average variation of $3$ mmHg in systole. The
reflections in the aortic flow predictions are also affected. Abdullateef et
al. [87] investigate the impact of fluctuating vessel radius on wave
reflections finding that as vessel radius decreases, wave reflections increase
and begin to overlap. Accurately placing nodes in tapering vessels is
therefore crucial as an incorrect placement can cause incorrect reflective
wave properties. Due to the well-known interplay between geometry and
hemodynamics, it is imperative to generate accurate vessel geometry in the
absence of ground truth data. Although this study uses the cylindrical
Navier$-$Stokes equations to model hemodynamics, considering the concept of
steady flow within a cylinder, namely Poiseuille’s law $q=\frac{\pi
pr^{4}}{8\mu l}$ [88, 89], allows us to gain insight into how radius and
length influence pressure and flow dynamics. Clearly, $q$ and $p$ are directly
impacted by $r$ and $l$ in Poiseuille’s law. Notably, $r$ is raised to the
fourth power, whereas $l$ is raised to the first power. Thus, we can expect
altering radius to have a greater influence on hemodynamics than length. Our
junction node modification process has a greater impact on a vessel’s length
than radius (Figure 3), leading to minimal differences in hemodynamics in some
vessels (Figure 8).

Change point analysis. We sequentially combine two methods for change point
analysis, utilizing R/changepoint to detect the optimal number of change
points [47] and R/segmented to place them [51]. Both methods are
computationally fast and their combination enables accurate and automated
assessment of vascular data. R/changepoint determines the optimal change point
count using mean, variance, and slope to detect changes in the data. However,
it is an approximation algorithm, only giving an estimated change point
placement rather than an exact location [47]. In contrast, R/segmented
accurately places change points through a customizeable and flexible
algorithm, where users define linear regression models based on prior
knowledge of data structure [51]. This package requires users to input the
number of change points to be detected within each dataset [51], so we use the
result from R/changepoint as an input. Coupling these two packages allows us
to create an automated pipeline for accurate vascular network analysis,
requiring little user input. This is advantageous when comparing multiple
vascular datasets, which may encompass thousands of vessels.

Limiting datasets to have up to three change points also provides efficient
and effective results. This ensures that detected change points capture
meaningful transitions in a vessel’s structure, such as junctions or image
resolution loss, while preventing overfitting of the data. It enables
consistent and computationally fast analysis, allowing for the development of
a robust algorithm for extracting representative radii (Section 2.5). While we
could define exactly three change points in all datasets, this leads to
inaccurate change point placement in vessels with an optimal count of less
than three, as observed in $9\%$ of pulmonary vessels and $33\%$ of aortic
vessels.

Radius sampling. Figure 9 shows that vessel radius has a large impact on
hemodynamics, following the pattern we expect from Poiseuille’s law. Given the
predicted pressure and flow values for various radii, a small change in radius
causes a large change in hemodynamics. For example, when we perturb the LCC’s
radius following a distribution $\mathcal{N}(0.36,0.01)$, we observe a 70 mmHg
range in systolic pressure. This is also shown in a study by Colebank et al.
[11], where even slight changes in vessel dimensions influence pressure and
flow. Secomb [90] notes that blood flow can be widely modulated by small
changes in diameter and emphasizes the need to precisely identify vessel
dimensions to understand how blood is distributed throughout the body . We
note that pulmonary blood pressures fluctuate throughout the respiratory
cycle, evident when using right-heart catheterization for pulmonary artery
measurements [91]. When measuring the same subject over time, a slight shift
in pressures is expected [92]. Our pressure range obtained from sampling radii
fall within the anticipated range, depending on how much radius is changed. It
is common practice to segment cardiovascular networks to build patient-
specific models used as non-invasive methods to monitor and understand disease
[42, 15, 41, 93, 94]. Our findings demonstrate that accurate radius
measurements are necessary to accurately tune models to patient data,
especially when using patient-specific parameters, such as vessel stiffness.

Limitations and future developments. This paper outlines a semi-automatic
framework to extract and refine patient-specific geometries from medical
images. We use a linear regression model to define change points for both non-
tapering and tapering vessels. With these change points, we are able to
reliably avoid unclear junction regions and detect representative radii.
However, since we use an exponential function to model tapering radii, we can
reduce computational complexity by directly defining change points with an
exponential model in those cases [95].

Another limitation is that we sample radii for each vessel independently.
While this allows us to investigate how geometry impacts hemodynamics, it
overlooks the intimate connection between vessels in the system. When we
sample independently, it is possible for daughter radii to exceed those of
their parent, violating the morphometric condition we enforce in Section 2.5.
In the future, we will implement a covariance structure that considers the
link between parents and daughter vessels when sampling [96].

We generate a relatively small number of segmentations for analysis in this
study due to the time-consuming nature of segmenting complex structures. We
plan to generate more segmentations of pulmonary and aortic vessels to further
test our method. To do so, we will use automatic segmentation software, such
as [97], and compare this to the manual segmentations we have already
completed.

##  5 Conclusion

This study demonstrates the need to investigate uncertainty in the network
generation process to create patient-specific vascular network models. We find
notable variation in networks, especially in the smallest vessels, when
segmenting the same scans multiple times. In addition, we find that relatively
small modulations in length and radius of vessels have large effects when
propagated to CFD models predicting hemodynamics. To minimize this
uncertainty, we develop a robust pipeline to determine radius and length of
vessels using change points. Based on our results, we emphasize the
significance of quantifying uncertainty when generating networks for use in
medical applications. Current diagnostic procedures of many diseases require
invasive protocol, including right-heart catheterization, but CFD modeling can
mitigate the need for this by providing insight in silico. Despite the
potential clinical advances CFD provides, medical conclusions should not be
made without first assessing the level of uncertainty present in diagnostic
information, fully understanding limitations and potential errors. In the
absence of fixed geometric values of vessels due to scanning ambiguity,
researchers must navigate uncertain measurements through the development of
effective methodologies, such as the one we have outlined here, for
increasingly accurate analyses.

Software and code implementing the developed methodology can be found at
github.com/msolufse.

MAB and AMT designed the study, developed algorithms, performed hemodynamic
simulations, and drafted the manuscript. DG, AJ, YL, ES, IS, and ZT
contributed to algorithm development and image segmentation. JDW and CP
provided aortic vasculature data and advised its use. DH provided statistical
expertise and conceived the study. MSO conceived and coordinated the study.

We have no competing interests.

This work was supported by the National Science Foundation (grant numbers
DGE-2137100, DMS-2051010), the National Security Agency (grant number
H98230-23-1-0009), and the National Institute of Health (grant numbers NIH-
HLBI RO1HL147590-02, NIH-5 T34 GM 131947-5). Any opinions, findings, and
conclusions expressed in this material are those of the authors and do not
necessarily reflect the views of the NSF, NSA, or NIH.

We thank Dr. Dan Lior for performing 4D-MRI registration and data extraction.

## References

  * [1] Y Chassidim, U Vazana, O Prager, R Veksler, G Bar-Klein, K Schoknecht, M Fassler, S Lublinsky and I Shelef  “Analyzing the blood–brain barrier: The benefits of medical imaging in research and clinical practice”  In _Semin. Cell Dev. Biol._ 38, 2015, pp. 43–52  DOI: https://doi.org/10.1016/j.semcdb.2014.11.007
  * [2] K Doi  “Diagnostic imaging over the last 50 years: research and development in medical imaging science and technology”  In _Phys. Med. Biol._ 51.13, 2006, pp. R5–R27  DOI: https://doi.org/10.1088/0031-9155/51/13/r02
  * [3] V S Khoo, D P Dearnaley, D J Finnigan, A Padhani, S F Tanner and M O Leach  “Magnetic resonance imaging (MRI): considerations and applications in radiotherapy treatment planning”  In _Radiother. Oncol._ 42.1, 1997, pp. 1–15  DOI: https://doi.org/10.1016/s0167-8140(96)01866-x
  * [4] H Alasti, Y B Cho, A D Vandermeer, A Abbas, B Norrlinger, S Shubbar and A Bezjak  “A novel four-dimensional radiotherapy method for lung cancer: imaging, treatment planning and delivery”  In _Phys. Med. Biol._ 51.12, 2006, pp. 3251–3267  DOI: https://doi.org/10.1088/0031-9155/51/12/017
  * [5] D L Pham, C Xu and J L Prince  “Current methods in medical image segmentation”  In _Annu. Rev. Biomed. Eng._ 2, 2000, pp. 315–337  DOI: https://doi.org/10.1146/annurev.bioeng.2.1.315
  * [6] D D Patil and S G Deore  “Medical image segmentation: a review”  In _Int. J. Comput. Sci. Mob. Computing_ 2, 2013, pp. 22–27
  * [7] N Sharma and L M Aggarwal  “Automated medical image segmentation techniques”  In _J. Med. Phys._ 35, 2010, pp. 3–14  DOI: https://doi.org/10.4103/0971-6203.58777
  * [8] L J O’Donnell  “Semi-automatic medical image segmentation”, 2001
  * [9] E M Van Rikxoort and B Van Ginneken  “Automated segmentation of pulmonary structures in thoracic computed tomography scans: a review”  In _Phys. Med. Biol._ 58, 2013, pp. R187–R220  DOI: https://doi.org/10.1088/0031-9155/58/17/r187
  * [10] T Buelow, R Wiemker, T Blaffert, C Lorenz and S Renisch  “Automatic extraction of the pulmonary artery tree from multi-slice CT data”  In _Medical Imaging 2005: Physiology, Function, and Structure from Medical Images_ 5746, 2005, pp. 730–740  SPIE  DOI: https://doi.org/10.1117/12.595286
  * [11] M J Colebank, M Paun, M U Quereshi, N C Chesler, D Husmeier, M S Olufsen and L Ellwein Fix  “Influence of image segmentation on one-dimensional fluid dynamics predictions in the mouse pulmonary arteries”  In _J. R. Soc. Interface_ 16, 2019, pp. 20190284  DOI: https://doi.org/10.1098/rsif.2019.0284
  * [12] P Lakhani  “The importance of image resolution in building deep learning models for medical imaging”  In _Radiol.: Artif._ 2.1 Radiological Society of North America, 2020, pp. e190177  DOI: https://doi.org/10.1148/ryai.2019190177
  * [13] R Mittal, J H Seo, V Vedula, Y J Choi, H Liu, H H Huang, S Jain, L Younes, T Abraham and R T George  “Computational modeling of cardiac hemodynamics: Current status and future outlook”  In _J. Comput. Phys._ 305 Elsevier, 2016, pp. 1065–1082  DOI: https://doi.org/10.1016/j.jcp.2015.11.022
  * [14] K M Chinnaiyan et al.  “Rationale, design and goals of the HeartFlow assessing diagnostic value of non-invasive FFRCT in Coronary Care (ADVANCE) registry”  In _J. Cardiovasc. Comput. Tomogr._ 11.1 Elsevier, 2017, pp. 62–67  DOI: https://doi.org/10.1016/j.jcct.2016.12.002
  * [15] A M Taylor-LaPole, M J Colebank, J D Weigand, M S Olufsen and C Puelz  “A computational study of aortic reconstruction in single ventricle patients”  In _Biomech. Model. Mechanobiol._ 22, 2023, pp. 357–377  DOI: https://doi.org/10.1007/s10237-022-01650-w
  * [16] M A Bartolo, M U Qureshi, M J Colebank, N C Chesler and M S Olufsen  “Numerical predictions of shear stress and cyclic stretch in pulmonary hypertension due to left heart failure”  In _Biomech. Model. Mechanobiol._ 21, 2022, pp. 363–381  DOI: https://doi.org/10.1007/s10237-021-01538-1
  * [17] C Battista, D Bia, YZ Germá, R L Armentano, M A Haider and M S Olufsen  “Wave propagation in a 1d fluid dynamics model using pressure-area measurements from ovine arteries”  In _J. Mech. Med. Biol._ 16.02 World Scientific, 2016, pp. 1650007  DOI: https://doi.org/10.1142/s021951941650007x
  * [18] L Formaggia, D Lamponi, M Tuveri and A Veneziani  “Numerical modeling of 1D arterial networks coupled with a lumped parameters description of the heart”  In _Comput. Methods Biomech. Biomed. Eng._ 9.5 Taylor & Francis, 2006, pp. 273–288  DOI: https://doi.org/10.1080/10255840600857767
  * [19] R A Gray and P Pathmanathan  “Patient-specific cardiovascular computational modeling: diversity of personalization and challenges”  In _J. Cardiovasc. Transl. Res._ 11 Springer, 2018, pp. 80–88  DOI: https://doi.org/10.1007/s12265-018-9792-2
  * [20] W Yang, M Dong, M Rabinovitch, F P Chan, A L Marsden and J A Feinstein  “Evolution of hemodynamic forces in the pulmonary tree with progressively worsening pulmonary arterial hypertension in pediatric patients”  In _Biomech. Model. Mechanobiol._ 18 Springer, 2019, pp. 779–796  DOI: https://doi.org/10.1007/s10237-018-01114-0
  * [21] R Botnar, G Rappitsch, M B Scheidegger, D Liepsch, K Perktold and P Boesiger  “Hemodynamics in the carotid artery bifurcation:: a comparison between numerical simulations and in vitro MRI measurements”  In _J. Biomech._ 33.2 Elsevier, 2000, pp. 137–144  DOI: https://doi.org/10.1016/s0021-9290(99)00164-5
  * [22] Y Shi, P Lawford and R Hose  “Review of zero-D and 1-D models of blood flow in the cardiovascular system”  In _Biomed. Eng. Online_ 10 Springer, 2011, pp. 1–38  DOI: https://doi.org/10.1186/1475-925x-10-33
  * [23] S M Moore, K T Moorhead, J G Chase, T David and J Fink  “One-dimensional and three-dimensional models of cerebrovascular flow”  In _J. Biomech. Eng._ 127, 2005  DOI: https://doi.org/10.1115/1.1894350
  * [24] P Reymond, F Perren, F Lazeyras and N Stergiopulos  “Patient-specific mean pressure drop in the systemic arterial tree, a comparison between 1-D and 3-D models”  In _J. Biomech._ 45.15 Elsevier, 2012, pp. 2499–2505  DOI: https://doi.org/10.1016/j.jbiomech.2012.07.020
  * [25] P J Blanco, C A Bulant, L O Müller, G D M Talou, C G Bezerra, P A Lemos and R A Feijoö  “Comparison of 1D and 3D models for the estimation of fractional flow reserve”  In _Sci. Rep._ 8.1 Nature Publishing Group UK London, 2018, pp. 17275  DOI: https://doi.org/10.1038/s41598-018-35344-0
  * [26] R Izzo, D Steinman, S Manini and L Antiga  “The vascular modeling toolkit: a python library for the analysis of tubular structures in medical images”  In _J. Open Source Softw._ 3.25, 2018, pp. 745  DOI: https://doi.org/10.21105/joss.00745
  * [27] L Ellwein, D S Marks, R Q Migrino, W D Foley, S Sherman and J F LaDisa Jr  “Image-based quantification of 3D morphology for bifurcations in the left coronary artery: Application to stent design”  In _Catheter. Cardiovasc. Interv._ 87.7 Wiley Online Library, 2016, pp. 1244–1255  DOI: https://doi.org/10.1002/ccd.26247
  * [28] A Updegrove, N M Wilson, J Merkow, H Lan, A L Marsden and S C Shadden  “SimVascular: an open source pipeline for cardiovascular simulation”  In _Ann. Biomed._ 45 Springer, 2017, pp. 525–541  DOI: https://doi.org/10.1007/s10439-016-1762-8
  * [29] D Chen and L D Cohen  “Piecewise geodesics for vessel centerline extraction and boundary delineation with application to retina segmentation”  In _Scale Space and Variational Methods in Computer Vision_ 9087, 2015, pp. 270–281  DOI: https://doi.org/10.1007/978-3-319-18461-6_22
  * [30] M Couprie and G Bertrand  “Asymmetric parallel 3D thinning scheme and algorithms based on isthmuses”  In _Pattern Recognit. Lett._ 76 Elsevier, 2016, pp. 22–31  DOI: https://doi.org/10.1016/j.patrec.2015.03.014
  * [31] P Hernandez-Cerdan  “Biopolymer networks: image analysis, reconstruction and modeling”, 2018
  * [32] M J Chambers  “Morphological Properties of Pulmonary Arteries in Control and Hypertensive Mice”  North Carolina State University, 2022
  * [33] S Sankaran, L Grady and C A Taylor  “Impact of geometric uncertainty on hemodynamic simulations using machine learning”  In _Comput. Meth. Appl. Mech. Eng_ 297, 2015, pp. 167–190  DOI: https://doi.org/10.1016/j.cma.2015.08.014
  * [34] P Chen, A Quarteroni and G Rozza  “Simulation-based uncertainty quantification of human arterial network hemodynamics”  In _Int. J. Numer. Methods Biomed. Eng._ 29.6 Wiley Online Library, 2013, pp. 698–721  DOI: https://doi.org/10.1002/cnm.2554
  * [35] G Bertaglia, A Navas-Montilla, A Valiani, M I M García, J Murillo and V Caleffi  “Computational hemodynamics in arteries with the one-dimensional augmented fluid-structure interaction system: viscoelastic parameters estimation and comparison with in-vivo data”  In _J. Biomech._ 100 Elsevier, 2020, pp. 109595  DOI: https://doi.org/10.1016/j.jbiomech.2019.109595
  * [36] J Gounley, M Vardhan and A Randles  “A computational framework to assess the influence of changes in vascular geometry on blood flow”  In _Proceedings of the Platform for Advanced Scientific Computing Conference_ , 2017, pp. 1–8  DOI: https://doi.org/10.1145/3093172.3093227
  * [37] M N Wilson, A K Ortiz and A B Johnson  “The vascular model repository: a public resource of medical imaging data and blood flow simulation results”  In _J. Med._ 7, 2013, pp. 040923  DOI: https://doi.org/10.1115/1.4025983
  * [38] A Federov et al.  “3D Slicer as an image computing platform for the Quantitative Imagine Network”  In _Magn. Reson. Imaging_ 20, 2012, pp. 1323–1341  DOI: https://doi.org/10.1016/j.mri.2012.05.001
  * [39] R Kikinis, S D Pieper and K Vosburgh  “3D Slicer: a platform for a subject-specific image analysis, visualization, and clinical support”  In _Intraoperative Imaging and Image-Guided Therapy_ New York, NY: Springer, 2014  DOI: https://doi.org/10.1007/978-1-4614-7657-3_19
  * [40] U Ayachit  “The ParaView Guide: A Parallel Visualization Application”  Clifton Park, NY, USA: Kitware, Inc., 2015
  * [41] L Antiga, M Piccinelli, L Botti, B Ene-Iordache, A Remuzzi and D A Steinman  “An image-based modeling framework for patient-specific computational hemodynamics”  In _Med. Biol. Eng. Comput._ 46, 2008, pp. 1097–1112  DOI: https://doi.org/10.1007/s11517-008-0420-1
  * [42] M J Colebank, M U Qureshi and M S Olufsen  “Sensitivity analysis and uncertainty quantification of 1-D models of pulmonary hemodynamics in mice under control and hypertensive conditions”  In _Int. J. Numer. Method. Biomed. Eng._ 37, 2021, pp. e3242  DOI: https://doi.org/10.1002/cnm.3242
  * [43] M Miller, N Johnston, I Livengood, M Spinelli, R Sazdanovic and M S Olufsen  “A topological data analysis study on murine pulmonary arterial trees with pulmonary hypertension”  In _Math. Biosci._ Elsevier, 2023, pp. 109056  DOI: https://doi.org/10.1016/j.mbs.2023.109056
  * [44] M S Olufsen, C S Peskin, W Y Kim, E M Pedersen, A Nadim and J Larsen  “Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions”  In _Ann. Biomed._ 28, 2000, pp. 1281–1299  DOI: https://doi.org/10.1114/1.1326031
  * [45] M J Chambers, M J Colebank, Qureshi M U, R Clipp and M S Olufsen  “Structural and hemodynamic properties of murine pulmonary arterial networks under hypoxia induced pulmonary hypertension”  In _Proc. Inst. Mech. Eng. H._ 234, 2020, pp. 1312–1329  DOI: https://doi.org/10.1177/0954411920944110
  * [46] M I Townsley  “Structure and composition of pulmonary arteries, capillaries and veins”  In _Compr. Physiol._ 2 NIH Public Access, 2012, pp. 675–709  DOI: https://doi.org/10.1002/cphy.c100081
  * [47] R Killick and I A Eckley  “changepoint: An R Package for Changepoint Analysis”  In _J. Stat. Softw._ 58, 2014, pp. 1–19  DOI: https://doi.org/10.18637/jss.v058.i03
  * [48] I A Eckley, P Fearnhead and R Killick  “Analysis of changepoint models”  In _Bayesian Time Series Modelcustomizeable statistical change point algorithm[50, 51]s_ Cambridge, UK: Cambridge University Press, 2011
  * [49] R Killick, P Fearnhead and I A Eckley  “Optimal detection of changepoints with a linear computational cost”  In _JASA_ 107.500 Taylor & Francis, 2012, pp. 1590–1598
  * [50] V Muggeo  “Estimating Regression Models with Unknown Break-Points”  In _Stat. Med._ 22, 2003, pp. 3055–3071  DOI: https://doi/org/10.1002/sim.1545
  * [51] V M R Muggeo  “Segmented: an R package to fit regression models with broken-line relationships”  In _R News_ 8, 2008, pp. 20–25
  * [52] C G Caro, T J Pedley, R C Schroter and W A Seed  “The Mechanics of the Circulation”  Cambridge, UK: Oxford Press, 1978
  * [53] H Cui, X Liu and N Huang  “Pulmonary vessel segmentation based on orthogonal fused U-Net${++}$ of chest CT images”  In _Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part VI 22_ , 2019, pp. 293–300  Springer  DOI: https://doi.org/10.1007/978-3-030-32226-7_33
  * [54] D Lesage, E D Angelini, I Bloch and G Funka-Lea  “A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes”  In _Med. Image Anal._ 13.6 Elsevier, 2009, pp. 819–845  DOI: https://doi.org/10.1016/j.media.2009.07.011
  * [55] M Zamir  “Nonsymmetrical bifurcations in arterial branching”  In _J. Gen. Physiol._ 72.6, 1978, pp. 837–845  DOI: https://doi.org/10.1085/jgp.72.6.837
  * [56] C D Murray  “The physiological principle of minimum work applied to the angle of branching of arteries”  In _J. Gen. Physiol._ 9.6 The Rockefeller University Press, 1926, pp. 835  DOI: https://doi.org/10.1085/jgp.9.6.835
  * [57] H B M Uylings  “Optimization of diameters and bifurcation angles in lung and vascular tree structures”  In _Bull. Math. Biol._ 39.5 Elsevier, 1977, pp. 509–520  DOI: https://doi.org/10.1016/s0092-8240(77)80054-2
  * [58] MathWorks  “MATLAB version: 9.13.0 (R2022b)”  Natick, Massachusetts, United States: The MathWorks Inc., 2022  URL: https://www.mathworks.com
  * [59] W W Nichols, M F O’Rourke and W L Kenney  “McDonald’s Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles”  Boca Raton, FL: CRC Press, 2011
  * [60] E M Pedersen, H W Sung, A C Burlson and A P Yoganathan  “Two-dimensional velocity measurements in a pulsatile flow model of the normal abdominal aorta simulating different hemodynamic conditions”  In _J. Biomech._ 26.10 Elsevier, 1993, pp. 1237–1247  DOI: https://doi.org/10.1016/0021-9290(93)90071-l
  * [61] A R Pries, D Neuhaus and P Gaehtgens  “Blood viscosity in tube flow: dependence on diameter and hematocrit”  In _Am. J. Physiol._ 263, 1992, pp. H1770–H1778  DOI: https://doi.org/10.1152/ajpheart.1992.263.6.h1770
  * [62] M U Qureshi, M J Colebank, L M Paun, L Ellwein Fix, N Chesler, M A Haider, N A Hill, D Husmeier and M S Olufsen  “Hemodynamic assessment of pulmonary hypertension in mice: a model-based analysis of the disease mechanism”  In _Biomech. Model. Mechanobiol._ 18 Springer, 2019, pp. 219–243  DOI: https://doi.org/10.1007/s10237-018-1078-8
  * [63] L M Paun, M J Colebank, M S Olufsen, N A Hill and D Husmeier  “Assessing model mismatch and model selection in a Bayesian uncertainty quantification analysis of a fluid-dynamics model of pulmonary blood circulation”  In _J. R. Soc. Interface_ 17.173, 2020, pp. 20200886  DOI: https://doi.org/10.1098/rsif.2020.0886
  * [64] J Kauffman and A G Schering  “Analysis of variance ANOVA”  In _WJSM_ , 2007  DOI: https://doi.org/10.1002/9780471462422.eoct017
  * [65] R Core Team  “R: A Language and Environment for Statistical Computing”, 2013  R Foundation for Statistical Computing  URL: http://www.R-project.org/
  * [66] Q R Grafton, H W Nelson, N R Lambie and P R Wyrwoll  “Coefficient of Variation”  In _The SAGE Encyclopedia of Research Design_ , 2022  URL: https://api.semanticscholar.org/CorpusID:199347749
  * [67] S Weglarczyk  “Kernel density estimation and its application”  In _ITM Web of Conferences_ 23, 2018, pp. 37
  * [68] A Norouzpour, Z Hooshyar and A Mehdizadeh  “Autoregulation of blood flow: vessel diameter changes in response to different temperatures”  In _J. Biomed. Phys. Eng._ 3.2 Shiraz University of Medical Sciences, 2013, pp. 63–66  URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4204495/
  * [69] J Vernon Jeffords and M H Knisely  “Concerning the geometric shapes of arteries and arterioles: A contribution to the biophysics of health, disease, and death”  In _Angiology_ 7.2 Sage Publications Sage CA: Thousand Oaks, CA, 1956, pp. 105–136  DOI: https://doi.org/10.1177/000331975600700202
  * [70] R E Feaver, B D Gelfand and B R Blackman  “Human haemodynamic frequency harmonics regulate the inflammatory phenotype of vascular endothelial cells”  In _Nat. Commun._ 4.1 Nature Publishing Group UK London, 2013, pp. 1525  DOI: https://doi.org/10.1038/ncomms2530
  * [71] E A V Jones, F Noble and A Eichmann  “What determines blood vessel structure? Genetic prespecification vs. hemodynamics”  In _Physiol. J._ 21.6 American Physiological Society, 2006, pp. 388–395  DOI: https://doi.org/10.1152/physiol.00020.2006
  * [72] S Lehoux and A Tedgui  “Cellular mechanics and gene expression in blood vessels”  In _J. Biomech._ 36.5 Elsevier, 2003, pp. 631–643  DOI: https://doi.org/10.1016/s0021-9290(02)00441-4
  * [73] D DeMers and D Wachs  “Physiology: Mean arterial pressure”  In _StatPearls Publishing_ 30, 2019, pp. 2019  URL: https://europepmc.org/article/nbk/nbk538226#_article-24804_s1_
  * [74] M Lins, J Vandevenne, M Thillai, B R Lavon, M Lanclus, S Bonte, R Godon, I Kendall, J De Backer and W De Backer  “Assessment of small pulmonary blood vessels in COVID-19 patients using HRCT”  In _Acad. Radiol._ 27.10 Elsevier, 2020, pp. 1449–1455  DOI: https://doi.org/10.1016/j.acra.2020.07.019
  * [75] S Hirsch, J Reichold, M Schneider, G Székely and B Weber  “Topology and hemodynamics of the cortical cerebrovascular system”  In _J. Cereb. Blood Flow Metab._ 32.6, 2012, pp. 952–967  DOI: https://doi.org/10.1038/jcbfm.2012.39
  * [76] P Horssen, J T B Lier, J P Wijngaard, E VanBavel, I E Hoefer, J A Spaan and Siebes M  “Influence of segmented vessel size due to limited imaging resolution on cornoary hyperemic flow prediction from arterial crown volume”  In _Am. J. Physiol._ 310, 2016, pp. H839–H846  DOI: https://doi.org/10.1152/ajpheart.00728.2015
  * [77] J C V Schwarz, M G Lier, J P H M Wijngaard, M Siebes and E VanBavel  “Topologic and hemodynamic characteristics of the human coronary arterial circulation”  In _Front. Physiol._ 10, 2020, pp. 1611  DOI: https://doi.org/10.3389/fphys.2019.01611
  * [78] C F Baumgartner, K C Tezcan, K Chaitanya, A M Hotker, U J Muehlematter, K Schawkat, A S Becker, O Donati and E Konukoglu  “Medical Image Computing and Computer Assisted Intervention”  Springer, 2019
  * [79] S Wang, B Peplinski, L Lu, W Zhang, J Liu, Z Wei and R M Summers  “Sequential Monte Carlo tracking for marginal artery segmentation on CT angiography by multiple cue fusion”  In _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2013: 16th International Conference, Nagoya, Japan, September 22-26, 2013, Proceedings, Part II 16_ , 2013, pp. 518–525  Springer  DOI: https://doi.org/10.1007/978-3-642-40763-5_64
  * [80] D Lesage, E D Angelini, G Funka-Lea and I Bloch  “Adaptive particle filtering for coronary artery segmentation from 3D CT angiograms”  In _Comput. Vis. Image Underst._ 151, 2016, pp. 29–46  DOI: https://doi.org/10.1016/j.cviu.2015.11.009
  * [81] A Mharib, A R Ramli, S Mashohor and R B Mahmood  “Survey on liver CT segmentation methods”  In _Artif. Intell. Rev._ 37, 2016, pp. 83–95  DOI: https://doi.org/10.1007/s10462-011-9220-3
  * [82] M M Hoeper et al.  “Definitions and diagnosis of pulmonary hypertension”  In _J. Am. Coll. Cardiol._ 62.25S, 2013, pp. D42–D50  DOI: https://doi.org/10.1016/j.jacc.2013.10.032
  * [83] L M Paun, A Borowska, M J Colebank, M S Olufsen and D Husmeier  “Inference in Cardiovascular Modelling Subject to Medical Interventions”  International Aset Inc., 2021  DOI: https://doi.org/10.11159/icsta21.109
  * [84] A Borowska, M J Colebank, M S Olufsen and D Husmeier  “Evaluating Closed-Loop Effects from Vasodilator Administration for Pulmonary Hypertension Treatment”  In _Proceedings of the 5th International Conference on Statistics: Theory and Applications_ , 2023  DOI: https://doi.org/10.11159/icsta23.162
  * [85] G E Ragkousis, N Curzen and N W Bressloff  “Multi-objective optimisation of stent dilation strategy in a patient-specific coronary artery via computational and surrogate modelling”  In _J. Biomech._ 49.2 Elsevier, 2016, pp. 205–215  DOI: https://doi.org/10.1016/j.jbiomech.2015.12.013
  * [86] S Epstein, M Willemet, P J Chowienczyk and J Alastruey  “Reducing the number of parameters in 1D arterial blood flow modeling: less is more for patient-specific simulations”  In _Am. J. Physiol. Heart Circ. Physiol._ 309.1 American Physiological Society Bethesda, MD, 2015, pp. H222–H234  DOI: https://doi.org/10.1152/ajpheart.00857.2014
  * [87] S Abdullateef, J Mariscal-Harana and A W Khir  “Impact of tapering of arterial vessels on blood pressure, pulse wave velocity, and wave intensity analysis using one-dimensional computational model”  In _Int. J. Numer. Meth. Biomed. Eng._ , 2020, pp. e3312  DOI: https://doi.org/10.1002/cnm.3312
  * [88] J Pfitzner  “Poiseuille and his law”  In _Anaesthesia_ 31, 1976, pp. 273–275  DOI: https://doi.org/10.1111/j.1365-2044.1976.tb11804.x
  * [89] D Demers and D Wachs  “StatPearls”  StatPearls Publishing, 2022
  * [90] T W Secomb  “Hemodynamics”  In _Compr. Physiol._ 6.2, 2016, pp. 975–1003  DOI: https://doi.org/10.1002/cphy.c150038
  * [91] G Kovacs, A Avian, M Pienn, R Naeije and H Olschewski  “Reading pulmonary vascular pressure tracings. How to handle the problems of zero leveling and respiratory swings”  In _Am. J. Respir. Crit. Care Med._ 190.3 American Thoracic Society, 2014, pp. 252–257  DOI: https://doi.org/10.1164/rccm.201402-0269pp
  * [92] S Magder  “Heart-Lung interaction in spontaneous breathing subjects: the basics”  In _Ann. Transl. Med._ 6.18 AME Publications, 2018, pp. 348  DOI: https://doi.org/10.21037/atm.2018.06.19
  * [93] J M Carson, C Roobottom, R Alcock and P Nithiarasu  “Computational instantaneous wave-free ratio (IFR) for patient-specific coronary artery stenoses using 1D network models”  In _Int. J. Numer. Meth. Biomed. Engng._ 35, 2019, pp. e3255  DOI: https://doi.org/10.1002/cnm.3255
  * [94] T Dobroserdova, S Simakov, T Gamilov, R Pryamonosov and E Sakharova  “Patient-specific blood flow modelling for medical applications”  In _MATEC Web of Conferences_ 76, 2016, pp. 05001  EDP Sciences  DOI: https://doi.org/10.1051/matecconf/20167605001
  * [95] J Chen and A K Gupta  “Parametric Statistical Change Point Analysis”  Boston, MA: Birkhauser, 2011
  * [96] A J McNeil  “Sampling nested Archimedean copulas”  In _J. Stat. Comput. Simul._ 78.6 TaylorFrancis, 2008, pp. 567–581  DOI: https://doi.org/10.1080/00949650701255834
  * [97] J Poletti, M Bach, S Yang, R Sexauer, B Stieltjes, D C Rotzinger, J Bremerich, A W Sauter and T Weikert  “Automated lung vessel segmentation reveals blood vessel volume redistribution in viral pneumonia”  In _Eur. J. Radiol._ 150 Elsevier, 2022, pp. 110259  DOI: https://doi.org/10.1016/j.ejrad.2022.110259