# Bayesian Windkessel calibration using optimized 0D surrogate models

Jakob Richter,1,2 Jonas Nitzler,2,3 Luca Pegolotti,1,4 Karthik Menon,1,4 Jonas
Biehler,2
Wolfgang A. Wall,2 Daniele E. Schiavazzi,5 Alison L. Marsden,1,4 Martin R.
Pfaller1,4

#### Affiliations

  1. 1.

Department of Pediatrics, Stanford University

  2. 2.

Institute for Computational Mechanics, Technical University of Munich

  3. 3.

Professorship for Data-driven Materials Modeling, Technical University of
Munich

  4. 4.

Institute for Computational and Mathematical Engineering, Stanford University

  5. 5.

Department of Applied and Computational Mathematics and Statistics, University
of Notre Dame

#### Abstract

Boundary condition (BC) calibration to assimilate clinical measurements is an
essential step in any subject-specific simulation of cardiovascular fluid
dynamics. Bayesian calibration approaches have been successfully employed for
quantifying the uncertainties inherent in identified parameters. Yet,
routinely estimating the posterior distribution for all BC parameters in 3D
simulations has been unattainable due to the infeasible computational demand.
We propose an efficient method to identify Windkessel parameter posteriors
using results from a single high-fidelity three-dimensional (3D) model
evaluation. We only evaluate the 3D model once for an initial choice of BCs
and use the result to create a highly accurate zero-dimensional (0D)
surrogate. We then perform Sequential Monte Carlo (SMC) using the optimized 0D
model to derive the high-dimensional Windkessel BC posterior distribution. We
validate this approach in a publicly available dataset of $N=72$ subject-
specific vascular models. We found that optimizing 0D models to match 3D data
a priori lowered their median approximation error by nearly one order of
magnitude. In a subset of models, we confirm that the optimized 0D models
still generalize to a wide range of BCs. Finally, we present the high-
dimensional Windkessel parameter posterior for different measured signal-to-
noise ratios in a vascular model using SMC. We further validate that the
0D-derived posterior is a good approximation of the 3D posterior. The minimal
computational demand of our method using a single 3D simulation, combined with
the open-source nature of all software and data used in this work, will
increase access and efficiency of Bayesian Windkessel calibration in
cardiovascular fluid dynamics simulations.

##  1 Introduction

Subject-specific Computational Fluid Dynamics (CFD) is an emerging tool in the
treatment of cardiovascular diseases [1]. Such CFD models require solving an
initial boundary value problem where the domain of interest is commonly
modeled with high-fidelity three-dimensional (3D) models while zero-
dimensional (0D) boundary conditions (BCs) physiologically approximate the
downstream vasculature (Figure 1, left). A common choice for BCs is the three-
element Windkessel model [2], which approximates the vascular resistance and
compliance of distal vascular beds that are not resolved in image data and
computationally intractable in 3D simulations. 0D models (Figure 1, middle and
right) can estimate bulk blood flow and pressure at a fraction of the
computational cost [3] of 3D simulations. In addition, 0D models are used as
efficient (albeit sometimes less accurate) substitutes for 3D models in the
simulation domain, e.g., to predict the pressure drop in aortic coarctations
[4] or coronary arteries [5]. They are especially useful in combination with
3D models in many-query settings, such as uncertainty quantification or
optimization problems involving thousands of model evaluations [6].

3D

model $\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})$

Geometric 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\boldsymbol{\alpha}_{\text{b}})$

Optimized 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$

Figure 1:  Illustration of different blood flow models used in this work.
Inflow (Inflow) and Windkessel boundary conditions (WK) with parameters
$\boldsymbol{\theta}$ are identical in all models. The geometry model is
separated into branch (light) and junction (dark) regions. The geometric 0D
model (middle) consists of BloodVessel elements (BV) for branches with
parameters $\boldsymbol{\alpha}_{\text{b}}$ derived purely from the 3D
geometry. The optimized 0D model (right) includes additional
BloodVesselJunction elements (BVJ) for junctions and all parameters
$\hat{\boldsymbol{\alpha}}$ are optimized from 3D simulation results.

Windkessel BC parameters are patient-specific, meaning they must be calibrated
for each patient based on clinical measurements [7], typically from 4D flow
magnetic resonance imaging (MRI) or Doppler ultrasound flow measurements and
catheter or cuff pressure measurements. In the following, these clinical
measurements are referred to as observations. These observations are usually
corrupted by noise, for example, from the limited accuracy of a measurement
device or the unpredictable positioning of a sensor. BCs estimated based on
these measurements are thus more accurately represented by a probability
distribution over a range of possible values rather than discrete parameter
values. The benefit of quantifying these inherent uncertainties is twofold.
First, it provides a measure of confidence when selecting a discrete point
estimate of BCs for further analysis. Second, we can further propagate the
uncertainties in BCs to uncertainties in model predictions for clinical
decision-making. In this work, we infer the posterior distribution of BC
parameters given noisy clinical observations. Our algorithm requires only
minimal 3D simulation costs, i.e., only a single fully resolved 3D simulation
over one cardiac cycle. To achieve this, we propose two novel techniques.
First, we present a fast deterministic optimization method to yield a highly
accurate 0D model from 3D data. Second, we use the 0D model to predict the
sought posterior distribution of the BC parameters given noisy observations.

The accuracy of 0D models compared to 3D fluid dynamics has been extensively
quantified in previous studies [5, 8, 9, 10, 3]. Automatic or semi-automatic
tools generate these 0D models from 3D geometries, which typically achieve
maximum approximation errors of bulk flow and pressure at their outlets of
around 10%, depending on the anatomy and the error metric used. In some
examples, resistances in 0D models have been calibrated semi-automatically or
automatically to improve agreement between 0D and 3D simulations. For
pulmonary models, Lan et al. [11] used a derivative-free Nelder-Mead
optimization, whereas Lee et al. [12] used iterative linear corrections. Tran
et al. [13] iteratively refined scaling global factors for key coronary model
outputs. In this work, we propose a deterministic least-squares procedure to
optimize a 0D model to match an initial 3D model. Our algorithm optimizes all
0D parameters in the model, i.e., resistances, inductances, capacitances, and
non-linear stenosis factors. Our 0D optimization method approximates the given
data in a least-squares sense, which makes it applicable to arbitrary vascular
anatomies. In this work, we demonstrate its use in aorta, aortafemoral,
coronary, and pulmonary models.

We refer to Marsden [14] for a general overview of optimization in
cardiovascular modeling and to Ninos et al. [15] for an overview of
uncertainty quantification (UQ) in hemodynamics. Furthermore, we point out
efficient developments in multi-fidelity uncertainty quantification [16, 17,
18, 19], which are especially useful for large-scale coupled problems, and can
drastically accelerate the otherwise expensive calculations. Different
strategies have been demonstrated to calibrate the parameters of Windkessel
BCs [20]. Spilker et al. [21] estimated Windkessel parameters by solving a
system of non-linear equations using a quasi-Newton method informed by a 0D
surrogate to compute the initial Jacobian. Ismail et al. [22] proposed an
adjoint-based approach to compute the initial Jacobian, thereby reducing the
sensitivity to accurate initial guesses. Tran et al. [13] and Menon et al.
[23, 24] used derivative-free optimization to estimate high-dimensional
Windkessel as well as other 0D BC parameters for coronary hemodynamics models
informed by clinical measurements of cardiac function and coronary perfusion.
However, deterministic optimization is prone to get stuck in unfavorable local
optima, especially in the presence of experimental noise [20]. More robust
approaches are based on Bayesian calibration, which interprets the calibration
as a statistical inference task, treating all variables as distributions.
Here, the solution is captured by a probability distribution for the
parameters. Pant et al. [25] demonstrated the application of Kalman filters
for Windkessel calibration using an unscented Kalman filter (UKF). The UKF was
applied to a 0D surrogate of the 3D model, which itself was subsequently
recalibrated based on 3D simulation results. Kalman filters sequentially
incorporate observations about a dynamic process (like time-dependent clinical
measurements) by propagating the mean and variance of the current parameter
estimate and correcting the latter based on the observations. The UKF is a
special form of Kalman filter, where dedicated control points are chosen for
the forward propagation to reduce the number of required model evaluations
[25]. Nevertheless, Kalman filters can only approximate Gaussian
distributions. It has been shown, however, that posterior distributions in
cardiovascular BCs deviate from Gaussian distributions, sometimes even
exhibiting multiple modes [26]. Other methods use multi-fidelity approaches
that rely on iteratively sampling the parameter space and learning a
relationship between high-fidelity 3D models and low-fidelity 0D or one-
dimensional (1D) models [27, 28].

In this work, we propose a new framework for calibrating Windkessel parameters
in cardiovascular models using Sequential Monte Carlo (SMC). As in Kalman
filter approaches, the calibration problem is stated as a Bayesian inference
task, and the result is given as the full posterior distribution over BC
parameters given the observations. The probabilistic approach also
incorporates prior knowledge about the BC parameters and information about the
noise structure in, e.g., clinical measurements used for the calibration. SMC
is employed to solve the Bayesian calibration problem. SMC can capture
arbitrary posteriors, including multimodal distributions. By exploring
intermediate distributions between the prior and the posterior, sampling is
sequentially refined in regions with high posterior density, leading to
computational savings. In our work, the workload of the Bayesian calibration
is completely shifted to 0D models. The 0D models are deterministically
optimized to a single 3D simulation, enabling the Bayesian calibration’s
accuracy. This yields highly accurate 0D models in unprecedented agreement
with the corresponding high-fidelity 3D models.

The remainder of this work is structured as follows: In Section 2, we outline
the basic concepts for 3D and 0D modeling and the optimization and Bayesian
calibration methods employed in this work. Section 3 presents the results of
the 0D model optimization and BC calibration algorithms and validates the
results against the high-fidelity 3D model. Section 4 summarizes the results
of this work and provides an outlook for future work.

##  2 Methods

In this section, we introduce Bayesian Windkessel calibration (2.1). Further,
we briefly review governing equations for high-fidelity three-dimensional (3D)
(2.2) and surrogate zero-dimensional (0D) (2.3) blood flow models. We then
extend these to the inverse problem of optimizing 0D models (2.4).

###  2.1 Bayesian Windkessel calibration with optimized 0D models

The goal of the Bayesian boundary condition (BC) calibration is to infer a
posterior distribution
$\mathfrak{p}_{\text{0D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}},\boldsymbol{\alpha})$
for the Windkessel BC parameters $\boldsymbol{\theta}$, informed by uncertain
clinical observations $\boldsymbol{y}_{\text{obs}}$. We perform this using a
0D model $\mathfrak{M}_{\text{0D}}$ with internal 0D element parameters
$\boldsymbol{\alpha}$. The distinction between the Windkessel BC parameters
$\boldsymbol{\theta}$ and the 0D element parameters $\boldsymbol{\alpha}$ is
crucial for the following discussion and will be explained in detail below. In
the following, we introduce probabilistic calibration (2.1.1), describe the
likelihood model selection in this work (2.1.2), and present the Sequential
Monte Carlo (SMC) algorithm (2.1.3).

####  2.1.1 Bayesian calibration

The posterior distribution reflects the solution of the calibration and can be
derived using Bayes’ rule [29]

$\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})=\frac{\mathfrak{p}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta})\cdot\mathfrak{p}(\boldsymbol{\theta})}{\mathfrak{p}(\boldsymbol{y}_{\text{obs}})}.$
(1)

The prior distribution $\mathfrak{p}(\boldsymbol{\theta})$ encodes the
knowledge about the parameters before new data is observed. The likelihood
$\mathfrak{p}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta})$ defines a
measure of proximity between the model response
$\boldsymbol{y}=\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})$ and the
observations $\boldsymbol{y}_{\text{obs}}$. The evidence
$\mathfrak{p}(\boldsymbol{y}_{\text{obs}})$ normalizes the former two
components to an actual probability density. This normalization constant is
typically difficult to evaluate, requiring the potentially high-dimensional
integration over $\boldsymbol{\theta}$. To avoid this complexity, many
algorithms make use of the unnormalized posterior

$\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})\propto\mathfrak{p}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta})\cdot\mathfrak{p}(\boldsymbol{\theta}),$
(2)

and the logarithmic version of the latter and conduct the normalization a
posteriori if needed.

Observations $\boldsymbol{y}_{\text{obs}}$ Likelihood
$\mathfrak{p}_{\text{3D}}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta})$
SMC iterations on $\boldsymbol{\theta}$ using $\mathfrak{M}_{\text{3D}}$ 3D
Model $\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})$ Prior
$\mathfrak{p}(\boldsymbol{\theta})$ Posterior
$\mathfrak{p}_{\text{3D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})$
Computationally infeasible Observations $\boldsymbol{y}_{\text{obs}}$
Likelihood
$\mathfrak{p}_{\text{0D}}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta},\boldsymbol{\alpha})$
SMC iterations on $\boldsymbol{\theta}$ using $\mathfrak{M}_{\text{0D}}$ 0D
Model $\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\boldsymbol{\alpha})$ Run
1) Geometric Prior $\mathfrak{p}(\boldsymbol{\theta})$ Posterior
$\mathfrak{p}_{\text{0D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}},\boldsymbol{\alpha})$
3D Model $\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{MAP}})$ Run 2)
LM optimization Point estimate $\boldsymbol{\theta}_{\text{MAP}}$ One-time
evaluation$\boldsymbol{\alpha}=\hat{\boldsymbol{\alpha}}$$\boldsymbol{\alpha}=\boldsymbol{\alpha}_{\text{b}}$Proposed
in this work$\approx$ Figure 2: Framework proposed in this work for Bayesian
Windkessel BC calibration, using 3D models $\mathfrak{M}_{\text{3D}}$
(expensive, red) and 0D models (cheap, green). The 0D models use geometric
parameters (gray) or optimized parameters (blue).

Bayesian calibration of Windkessel BC parameters $\boldsymbol{\theta}$ from
observations $\boldsymbol{y}_{\text{obs}}$ is a computationally challenging
task. Using SMC, it easily requires thousands of model evaluations
$\mathfrak{M}$, each requiring a 3D FEM solution (Figure 1, left), to yield
the posterior distribution
$\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})$ from a prior
distribution $\mathfrak{p}(\boldsymbol{\theta})$. Figure 2 (top) highlights
this process using the 3D model $\mathfrak{M}_{\text{3D}}$ (red) to obtain the
posterior distribution
$\mathfrak{p}_{\text{3D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})$.
While this would give an exact representation of the posterior distribution,
it is generally infeasible due to its excessive computational demand.

In this work, we combine models $\mathfrak{M}_{\text{0D}}$ and
$\mathfrak{M}_{\text{3D}}$ in a novel calibration method visualized in Figure
2 (bottom). Here, expensive statistical model evaluation is performed
exclusively with fast 0D models $\mathfrak{M}_{\text{0D}}$ (green), relying on
two crucial assumptions. First, we assume that we can construct highly
accurate 0D surrogate models that approximate
$\mathfrak{M}_{\text{3D}}\approx\mathfrak{M}_{\text{0D}}$. Second, we assume
that the 0D-derived posterior approximates
$\mathfrak{p}_{\text{3D}}\approx\mathfrak{p}_{\text{0D}}$. While 0D models are
orders of magnitude faster to evaluate, they are, in general, not sufficiently
accurate to justify these assumptions. As part of our calibration process, we
thus optimize the 0D element parameters $\boldsymbol{\alpha}$ from a single 3D
simulation to greatly enhance 0D accuracy. This also allows us to demonstrate
the validity of the above two assumptions.

We perform two SMC calibration runs in the proposed calibration workflow. The
first SMC run (gray) is performed using a geometric 0D model (Figure 1,
middle)
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\boldsymbol{\alpha}_{\text{b}})$
composed of BloodVessel elements with parameters
$\boldsymbol{\alpha}_{\text{b}}$ derived purely from the blood vessel’s
geometry, containing the vessel branches’ resistances, inductances,
capacitances, and non-linear stenosis factors (Section 2.3) [3]. This yields a
first estimate for the posterior distribution
$\mathfrak{p}_{\text{0D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}},\boldsymbol{\alpha}_{\text{b}})$.
Using the parameter set $\boldsymbol{\theta}_{\text{MAP}}$ with the highest
posterior density, the so-called maximum a posteriori (MAP) estimate, we
perform one evaluation of the 3D model
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{MAP}})$ (Section 2.2).
These results are then used to determine the optimized 0D model (Figure 1,
right) parameters $\hat{\boldsymbol{\alpha}}$, including BloodVesselJunction
elements, using the Levenberg-Marquardt (LM) algorithm (Section 2.4). This
optimization yields a highly accurate 0D model to approximate
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{MAP}})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{MAP}},\hat{\boldsymbol{\alpha}})$.
Finally, we perform a second SMC (blue) run with the optimized 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$ to
derive the final posterior distribution
$\mathfrak{p}_{\text{0D}}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$
for the Windkessel BC parameters $\boldsymbol{\theta}$. Throughout this
manuscript, we refer to obtaining the optimal point estimate 0D parameters
$\hat{\boldsymbol{\alpha}}$ using LM as optimization. This is in contrast to
the Bayesian calibration of Windkessel parameters $\boldsymbol{\theta}$ with
posterior distribution
$\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})$ under
uncertain observations $\boldsymbol{y}_{\text{obs}}$ using SMC.

####  2.1.2 Zero-dimensional likelihood model

The forward model can only be evaluated pointwise, i.e., for a given set of
parameters $\boldsymbol{\theta}$. The likelihood model scores the model output
vector $\boldsymbol{y}=\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})$ for a
given parameter vector $\boldsymbol{\theta}$ based on the output’s proximity
to the observations $\boldsymbol{y}_{\text{obs}}$. Moreover, it incorporates
knowledge about the noise in the observations resulting from the clinical
measurements [30]. In this work, an additive noise model for the observations
is chosen according to Kaipio et al. [30] as

$\boldsymbol{y}_{\text{obs}}=\boldsymbol{y}+\boldsymbol{\varepsilon}=\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})+\boldsymbol{\varepsilon}.$
(3)

The noise is assumed to be normally distributed with zero mean and covariance
matrix $\boldsymbol{\Sigma}$, i.e.,
$\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma})$.
The covariance matrix is further assumed to be diagonal wherein the variance
vector $\boldsymbol{\sigma}^{2}$ denotes the diagonal entries of
$\boldsymbol{\Sigma}$. We choose the commonly used Gaussian likelihood model
here, which is defined as:

$\mathfrak{p}\left(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta}\right)=\frac{1}{\sqrt{\left(2\pi\right)^{n}\mathrm{det}\,\boldsymbol{\Sigma}}}\exp\left[-\frac{1}{2}\left(\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})-\boldsymbol{y}_{\text{obs}}\right)^{\top}\boldsymbol{\Sigma}^{-1}\left(\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})-\boldsymbol{y}_{\text{obs}}\right)\right]$
(4)

It should be noted that
$\mathfrak{p}\left(\boldsymbol{y}_{\text{obs}}|\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta})\right)$
is not a valid probability density in $\boldsymbol{\theta}$. Therefore, it is
often referred to as likelihood function in $\boldsymbol{\theta}$.

####  2.1.3 Sequential Monte Carlo

Different methods have been developed to solve the Bayesian calibration
problem. We employ the SMC algorithm, also known as particle filters [31, 32],
which is described in more detail in Appendix C. The posterior distribution
$\pi$ in SMC is approximated by weighted particles according to

$\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})\approx\sum_{i=1}^{k}W_{n+1}^{(i)}\,\delta_{\boldsymbol{\theta}_{n+1}^{(i)}}(\boldsymbol{\theta}),$
(5)

with normalized weights $W$, number of particles $k$, and final iteration
number $n$. Each particle $\delta_{\boldsymbol{\theta}_{n}^{(i)}}$ refers to a
Dirac delta distribution at particle location $\boldsymbol{\theta}_{i}$. Using
the parameter vector $\boldsymbol{\theta}$ with the highest posterior
probability [30], the so-called MAP estimate can be computed as

$\boldsymbol{\theta}_{\text{MAP}}=\mathrm{arg}\max_{\boldsymbol{\theta}}\left[\,\mathfrak{p}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})\,\right].$
(6)

In SMC, this is equivalent to selecting the particle location
$\boldsymbol{\theta}_{i}$ of the particle with the highest weight.

###  2.2 Three-dimensional blood flow modeling

The 3D model $\mathfrak{M}_{\text{3D}}$ models blood as a Newtonian fluid in
the rigid fluid domain $\Omega^{\text{3D}}$ over time $T$ (typically several
cardiac cycles) governed by the incompressible Navier-Stokes equations,

$\displaystyle\nabla\cdot\boldsymbol{\sigma}(\mathbf{u},p)$
$\displaystyle=\rho\left[\dot{\mathbf{u}}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right],\quad$
$\displaystyle\mathbf{x}\in\Omega^{\text{3D}},\quad t\in[0,T],$ (7)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,\quad$
$\displaystyle\mathbf{x}\in\Omega^{\text{3D}},\quad t\in[0,T],$ (8)

with velocity $\mathbf{u}$, pressure $p$, Cauchy stress tensor
$\boldsymbol{\sigma}=\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}})-p\mathbf{I}$,
constant density $\rho$, and constant dynamic viscosity $\mu$. Throughout this
manuscript, $\dot{\phantom{m}}$ denotes a time derivative. We prescribe BCs

$\displaystyle\mathbf{u}(\mathbf{x},t)$ $\displaystyle=\mathbf{0},\quad$
$\displaystyle\mathbf{x}\in\Gamma_{\text{wall}},\quad$ $\displaystyle
t\in[0,T],$ (9) $\displaystyle\mathbf{u}(\mathbf{x},t)$
$\displaystyle=\mathbf{u}_{\text{in}}(\mathbf{x},t),\quad$
$\displaystyle\mathbf{x}\in\Gamma_{\text{in}},\quad$ $\displaystyle
t\in[0,T],$ (10)
$\displaystyle\mathbf{f}_{\text{WK3}}\left(P_{i}(t),Q_{i}(t),\,\boldsymbol{\theta}_{i}\right)$
$\displaystyle=\mathbf{0},\quad$
$\displaystyle\mathbf{x}\in\Gamma_{\text{out},i},\quad$ $\displaystyle
t\in[0,T],$ (11) $\displaystyle\text{with\leavevmode\nobreak\
}p(\mathbf{x},t)$ $\displaystyle=P_{i}(t),$ (12) $\displaystyle Q_{i}(t)$
$\displaystyle=\int_{\Gamma_{\text{out},i}}\mathbf{u}(\mathbf{x},t)$
$\displaystyle\cdot\mathbf{n}\,\text{d}\mathbf{x},$ (13)

with a no-slip BC at the vessel wall $\Gamma_{\text{wall}}$, time-dependent
parabolic inflow profile $\mathbf{u}_{\text{in}}$ at the inlet
$\Gamma_{\text{in}}$. At each outlet $\Gamma_{\text{out},i}$, we prescribe a
three-element Windessel BC with parameters $\boldsymbol{\theta}_{i}$ that
relate flow rate $Q_{i}$ and the pressure $P_{i}$ at outlet $i$ through an
ordinary differential equation $\mathbf{f}_{\text{WK3}}$. Throughout this
work, we use the three-element Windkessel model defined in Appendix B.5 for
all outlets. We note that this choice does not restrict the BC calibration
approach proposed in this work as it applies generally to any BC type. For
details on the 0D-3D coupling approach, we refer to Vignon-Clementel et al.
[33, 34, 35]. The initial conditions are

$\displaystyle\mathbf{u}(\mathbf{x},t=0)$
$\displaystyle=\mathbf{u}_{0}(\mathbf{x}),\qquad$
$\displaystyle\mathbf{x}\in\Omega^{\text{3D}}$ (14) $\displaystyle
p(\mathbf{x},t=0)$ $\displaystyle=p_{0}(\mathbf{x}),\qquad$
$\displaystyle\mathbf{x}\in\Omega^{\text{3D}}$ (15)

with initial velocity field $\mathbf{u}_{0}$ and initial pressure field
$p_{0}$, which we estimate by projecting a periodic 0D solution to the 3D
geometry [36].

###  2.3 Zero-dimensional blood flow modeling

Analogous to electrical circuits, the elemental building blocks for blood
vessels of the 0D model $\mathfrak{M}_{\text{0D}}$ are defined as

$\displaystyle\Delta P=\mathcal{R}Q,\quad\Delta P=\mathcal{L}\dot{Q},\quad
Q=\mathcal{C}\dot{P},\quad\Delta P=\mathcal{S}Q|Q|,$ (16)

with 0D parameters resistance $\mathcal{R}$, inductance $\mathcal{L}$,
capacitance $\mathcal{C}$, and stenosis coefficient $\mathcal{S}$. Note that
only the stenosis term depends non-linearly on the flow $Q$. The operator,
$\Delta\bullet$ characterizes the difference between element outlet
$\bullet_{\text{out}}$ and inlet $\bullet_{\text{in}}$.

We use the method outlined in Pfaller et al. [3] to split the vessel into
branches and junctions. We combine the $\mathcal{RLCS}$ elements in (16) into
a BloodVessel lumped-parameter network (LPN) element as described in Appendix
B.3. We predict the branch parameter values $\boldsymbol{\alpha}_{\text{b}}$
of each BloodVessel element using information from the blood vessel geometry,
such as vessel length and cross-sectional area throughout the branch. In the
following sections, these models are referred to as geometric 0D models. We
also define a BloodVesselJunction LPN element to connect an inlet to two or
more outlets using the $\mathcal{RLS}$ elements with parameters
$\boldsymbol{\alpha}_{\text{j}}$, as described in Appendix B.4. The
BloodVesselJunction LPN element does not include a capacitance, which
facilitates the derivation of the Jacobian for the zero-dimensional model
optimization. Note that we do not currently infer BloodVesselJunction element
parameters from the blood vessel geometry. Predicting junction parameters from
the vessel geometry remains an active research topic [37, 38]. Thus, all
junction parameters $\boldsymbol{\alpha}_{\text{j}}$ are set to zero in
geometric 0D models, i.e., $\boldsymbol{\alpha}_{\text{j}}=\mathbf{0}$.
Nevertheless, we optimize parameters for both BloodVessel and
BloodVesselJunction elements in Section 2.4 using 3D simulation results. These
models are referred to as optimized 0D models.

Connecting all elements $e$ within a vascular tree forms a LPN with parameters
$\boldsymbol{\alpha}=\boldsymbol{\alpha}_{\text{b}}\cup\boldsymbol{\alpha}_{\text{j}}$.
As in the 3D model (Section 2.2), we prescribe the flow rate $Q$ at the inlet
and couple the outlets to Windkessel BCs. We assemble all elements $e$ and BCs
in our LPN into a global residual vector $\mathbf{r}$. The global system of
differential-algebraic equations can be split by the dependence on the
solution vector $\mathbf{y}$ and its time derivative $\dot{\mathbf{y}}$ as

$\displaystyle\mathbf{r}(\boldsymbol{\alpha},\mathbf{y},\dot{\mathbf{y}})=\mathbf{E}(\boldsymbol{\alpha})\cdot\dot{\mathbf{y}}+\mathbf{F}(\boldsymbol{\alpha})\cdot\mathbf{y}+\mathbf{c}(\boldsymbol{\alpha},\mathbf{y},\dot{\mathbf{y}}),$
(17)

with constant parameter-dependent system matrices $\mathbf{E}$ and
$\mathbf{F}$ and system vector $\mathbf{c}$ that depends non-linearly on
$\mathbf{y}$ and $\dot{\mathbf{y}}$. We compute pressure and flow rate
associated with a given set of 0D parameters $\boldsymbol{\alpha}$ by solving

$\displaystyle\text{Solve\leavevmode\nobreak\
}\mathbf{r}(\boldsymbol{\alpha},\mathbf{y},\dot{\mathbf{y}})\overset{!}{=}\mathbf{0}\text{\leavevmode\nobreak\
for\leavevmode\nobreak\ }\mathbf{y},\dot{\mathbf{y}}\text{\leavevmode\nobreak\
with given\leavevmode\nobreak\
}\boldsymbol{\alpha},\quad\mathbf{y}\in\Omega^{\text{0D}},\quad t\in[0,T],$
(18)

which, in the following, is referred to as the forward 0D problem. The node
points of the LPN, i.e., where 0D elements connect, make up the domain
$\Omega^{\text{0D}}$. We generate maps
$\Omega^{\text{0D}}\to\Omega^{\text{3D}}$ and
$\Omega^{\text{3D}}\to\Omega^{\text{0D}}$ through a region-growing algorithm
and centerline integration, respectively, to map results between 0D and 3D
models [36]. We prescribe BCs at the inlet and outlets

$\displaystyle Q_{j}(t)$ $\displaystyle=Q_{\text{in}}(t),\quad$ $\displaystyle
y_{j}\in\Gamma_{\text{in}},\quad$ $\displaystyle t\in[0,T],$ (19)
$\displaystyle\mathbf{f}_{\text{WK3}}\left(P_{j}(t),Q_{j}(t),\,\boldsymbol{\theta}_{j}\right)$
$\displaystyle=\mathbf{0},\quad$ $\displaystyle
y_{j}\in\Gamma_{\text{out},i},\quad$ $\displaystyle t\in[0,T],$ (20)

with inflow rate $Q_{\text{in}}$ and three-element Windkessel BCs
$\mathbf{f}_{\text{WK3}}$ corresponding to the 3D model inlet and outlet BCs
$\boldsymbol{\theta}$, (10) and (11), respectively. Furthermore, we prescribe
initial conditions

$\displaystyle\mathbf{y}(t=0)$ $\displaystyle=\bar{\mathbf{y}},\quad$
$\displaystyle\mathbf{y}\in\Omega^{\text{0D}},$ (21)

where we obtain $\bar{\mathbf{y}}$ from a steady-state simulation with mean
inflow and steady-state BCs. This lowers the number of cardiac cycles required
to reach a periodic solution compared to initializing with
$\mathbf{y}(t=0)=\mathbf{0}$. We use the implicit generalized-$\alpha$ method
to discretize the forward problem (18) in time and the Newton-Raphson method
to solve for the solution vector $\mathbf{y}$ in each time step. The details
of the method are given in Appendix B.1. Reusing terms from the residual (17),
the Jacobian matrix $\mathbf{K}$ for the time-discrete forward problem (18) is
computed as

$\displaystyle\mathbf{K}=\left.\frac{\partial\mathbf{r}}{\partial\dot{\mathbf{y}}}\right|_{\boldsymbol{\alpha}=\text{const}}=\alpha_{m}\left(\mathbf{E}+\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)+\alpha_{f}\gamma\Delta
t_{n}\left(\mathbf{F}+\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right),$
(22)

with time step $\Delta t_{n}$ and time integration parameters $\alpha_{f}$,
$\alpha_{m}$, and $\gamma$. All system matrices for the BloodVessel,
BloodVesselJunction, and Windkessel elements can be derived analytically and
are provided in Appendices B.3, B.4, and B.5, respectively.

###  2.4 Zero-dimensional model optimization

The proposed method for Bayesian Windkessel calibration relies on highly
accurate 0D models $\mathfrak{M}_{\text{0D}}$ that act as a surrogate for a
given 3D model $\mathfrak{M}_{\text{3D}}$. In general, 0D models with
geometric parameters $\boldsymbol{\alpha}_{\text{b}}$ are insufficiently
accurate. In this section, we describe how to derive an optimized 0D parameter
set $\hat{\boldsymbol{\alpha}}$ that better approximates
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{MAP}})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{MAP}},\hat{\boldsymbol{\alpha}})$
for a given set of Windkessel parameters $\boldsymbol{\theta}_{\text{MAP}}$.
Solving an inverse problem, we infer 0D model parameters
$\hat{\boldsymbol{\alpha}}$ that optimally (in the least-squares sense)
approximate given 0D solution vectors $\mathbf{y}$ and $\dot{\mathbf{y}}$. In
this work, $\mathbf{y}$ is extracted from a 3D model
$\mathfrak{M}_{\text{3D}}$ evaluated at the parameter set
$\boldsymbol{\theta}_{\text{MAP}}$. However, other sources like in vivo or in
vitro measurements are also possible. The time derivative $\dot{\mathbf{y}}$
can be either taken from $\mathfrak{M}_{\text{3D}}$ or approximated from
$\mathbf{y}$. We use the same residual $\mathbf{r}$ from (17) as in the
forward problem (18):

$\displaystyle\text{Solve\leavevmode\nobreak\
}\mathbf{r}(\boldsymbol{\alpha},\mathbf{y},\dot{\mathbf{y}})\overset{!}{=}\mathbf{0}\text{\leavevmode\nobreak\
for\leavevmode\nobreak\ }\boldsymbol{\alpha}\text{\leavevmode\nobreak\ with
given\leavevmode\nobreak\
}\mathbf{y},\dot{\mathbf{y}},\quad\mathbf{y}\in\Omega^{\text{0D}}.$ (23)

To solve the inverse problem, we formulate (23) as a least-squares problem

$\displaystyle\hat{\boldsymbol{\alpha}}=\arg\min_{\boldsymbol{\alpha}}S,\quad\text{with\leavevmode\nobreak\
}S=\sum_{j}^{D}r_{j}^{2}(\boldsymbol{\alpha},y_{j},\dot{y}_{j}),$ (24)

where $\hat{\boldsymbol{\alpha}}$ is the optimal parameter set that minimizes
the sum of squared residuals $S$, and $D$ is the number of observations. Here,
we solve for both branch and junction 0D parameters
$\hat{\boldsymbol{\alpha}}_{\text{b}}$ and
$\hat{\boldsymbol{\alpha}}_{\text{j}}$, respectively. The value $D$ is
typically given from the number of LPN node points times the number of
observed time steps. Note that, in general, any combination of temporal and
spatial observations is possible as long as the point estimate parameters
$\hat{\boldsymbol{\alpha}}$ remain deterministically identifiable. We use the
LM algorithm to solve the least squares problem (24) iteratively for the
optimal parameter set $\hat{\boldsymbol{\alpha}}$, which is described in more
detail in Appendix B.2. We calculate the Jacobian for each 0D element as the
gradient of the residual $\mathbf{r}$ with respect to the 0D parameters
$\boldsymbol{\alpha}$,

$\displaystyle\mathbf{J}=\left.\frac{\partial\mathbf{r}}{\partial\boldsymbol{\alpha}}\right|_{\mathbf{y},\dot{\mathbf{y}}=\text{const}}=\frac{\partial\mathbf{E}}{\partial\boldsymbol{\alpha}}\cdot\dot{\mathbf{y}}+\frac{\partial\mathbf{F}}{\partial\boldsymbol{\alpha}}\cdot\mathbf{y}+\frac{\partial\mathbf{c}}{\partial\boldsymbol{\alpha}}.$
(25)

Note that the solution vectors $\mathbf{y},\dot{\mathbf{y}}$ are constant in
(25), and the Jacobian can, again, be computed analytically. These terms are
straightforward to derive and are given for the BloodVessel and
BloodVesselJunction LPN elements in Appendices B.3 and B.4, respectively. The
resulting 0D models $\mathfrak{M}_{\text{0D}}$ with optimal parameters
$\hat{\boldsymbol{\alpha}}=\hat{\boldsymbol{\alpha}}_{\text{b}}\cup\hat{\boldsymbol{\alpha}}_{\text{j}}$
are denoted as optimized 0D models.

###  2.5 Open-source software

All results in this work were computed with open-source software. The 3D
simulations were computed in SimVascular’s finite element solver svSolver
[39]; see Esmaily Moghadam et al. [35, 40] for details. We further developed
svZeroDSolver [41], containing the algorithmic implementation of the methods
described in 2.3 and 2.4. The forward (18) and inverse problems (23) use the
same residual vector $\mathbf{r}$ (17) with linearizations with respect to the
solution vector $\mathbf{y}$ and the parameter vector $\boldsymbol{\alpha}$,
respectively. The forward problem requires time discretization, whereas the
inverse problem is simultaneously solved for all time steps. Nevertheless, the
shared residual allows for reuse of large parts of the code for both problems.
The system matrices for all 0D elements for forward and inverse problems are
given in Appendix B. The core functionality is implemented in C++ and based on
the fast linear algebra package Eigen. Python and C++ interfaces enable
seamless integration into other software [12, 24]. The calibration framework
svSuperEstimator [42] implements the Windkessel calibration described in
Section 2.1 and interfaces different tools for performing 0D and 3D
simulations, postprocessing routines and SMC Windkessel calibration. The SMC
implementation is based on the Python package particles [43]. All 72 vascular
models are publicly available from the Vascular Model Repository [44].

##  3 Computational results and discussion

Our results showcase the two-run process of Bayesian calibration of Windkessel
boundary condition (BC) parameters $\boldsymbol{\theta}$ given uncertain
observations $\boldsymbol{y}_{\text{obs}}$. First, Section 3.1 demonstrates
the approximation of three-dimensional (3D) models by zero-dimensional (0D)
models
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$
with optimized 0D parameters $\hat{\boldsymbol{\alpha}}$. Second, Section 3.2
provides posterior distributions for BC parameters $\boldsymbol{\theta}$
derived with Sequential Monte Carlo (SMC) using the optimized 0D models,
$\mathfrak{p}_{\text{3D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}})\approx\mathfrak{p}_{\text{0D}}(\boldsymbol{\theta}|\boldsymbol{y}_{\text{obs}},\hat{\boldsymbol{\alpha}})$.

###  3.1 Optimization of zero-dimensional models

In this section, we investigate the Levenberg-Marquardt (LM) 0D model
optimization performance with the algorithm presented in Section 2.4. The LM
optimizer derives new blood vessel and junction parameters
$\hat{\boldsymbol{\alpha}}$ from a 3D simulation result, enabling the
approximation
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$.
The 0D models are improved in two aspects: enhancing the geometrically derived
BloodVessel elements $\hat{\boldsymbol{\alpha}}_{\text{b}}$, and identifying
suitable parameters for BloodVesselJunction elements
$\hat{\boldsymbol{\alpha}}_{\text{j}}$ (which are zero in geometric 0D
models). We obtain the time derivative $\dot{\mathbf{y}}$ by fitting a
periodic cubic spline to $\mathbf{y}$, taking its derivative, and resampling
the solution to 100 time points. The number of data points $D$ was thus 100
times the number of lumped-parameter network (LPN) node points, which ranged
from 28 in aortas to 1460 in pulmonary arteries (mean 274).

To demonstrate the robustness of the 0D model optimization, we demonstrate the
results on a large dataset of 72 vascular models, each with given Windkessel
BC parameters $\boldsymbol{\theta}_{\text{VMR}}$ that were calibrated in prior
studies. This openly available dataset [44] is taken from Pfaller et al. [3]
and visualized in Figure 7. It includes diverse vascular anatomies, including
aortas, aortofemoral, pulmonary, and coronary arteries. Vascular states
include normal, (artificial) aortic coarctation, Kawasaki disease, single
ventricle disease, Marfan syndrome, abdominal aortic aneurysm, pulmonary
artery hypertension, and others. Note that the 3D dataset has rigid walls,
corresponding to a theoretical 0D capacitance of $\mathcal{C}=0$. However, we
found that using a small non-zero value for $\mathcal{C}$, thus modeling a
nearly rigid wall, is beneficial for 0D numerical performance [3]. For
completeness, we still report results for $\mathcal{C}$ in this section
despite it not having a physical meaning in our 0D models.

We first verify in Section 3.1.1 that we can identify 0D branch parameters
$\hat{\boldsymbol{\alpha}}_{\text{b}}$ from an optimization given ground truth
0D models with geometric branch parameters $\boldsymbol{\alpha}_{\text{b}}$.
We then demonstrate in Section 3.1.2 the accuracy of 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})$
with parameters $\hat{\boldsymbol{\alpha}}$ optimized from 3D models
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{VMR}})$. Finally, we
verify in Section 3.1.3 that the optimized 0D models generalize well to a wide
range of BC parameters $\boldsymbol{\theta}$. These results confirm that we
can confidently use optimized 0D models for the approximation
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$
in the Bayesian Windkessel calibration in Section 3.2.

####  3.1.1 Identifying zero-dimensional parameters from ground truth data

As a first step in assessing the 0D element optimization proposed in Section
2.4, we verify that our optimization can successfully recover ground truth 0D
geometric branch parameters $\boldsymbol{\alpha}_{\text{b}}$ from
$\mathfrak{M}_{\text{0D}}$ data:

$\displaystyle\text{Compare\leavevmode\nobreak\
}\hat{\boldsymbol{\alpha}}_{\text{b}}\text{\leavevmode\nobreak\
vs.\leavevmode\nobreak\
}\boldsymbol{\alpha}_{\text{b}}\text{\leavevmode\nobreak\
with\leavevmode\nobreak\
}\hat{\boldsymbol{\alpha}}_{\text{b}}=\arg\min_{\boldsymbol{\alpha}}\|\mathbf{r}(\boldsymbol{\alpha},\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\boldsymbol{\alpha}_{\text{b}}))\|^{2}.$
(26)

In the absence of any 3D-0D model discrepancy, the optimization should be able
to perfectly recover the ground truth 0D parameters
$\hat{\boldsymbol{\alpha}}_{\text{b}}\approx\boldsymbol{\alpha}_{\text{b}}$.
Recall from Section 1 that $\boldsymbol{\alpha}_{\text{j}}=\mathbf{0}$ in all
geometric models which results in optimized junction parameters
$\hat{\boldsymbol{\alpha}}_{\text{j}}=\mathbf{0}$. We ran forward models
$\mathfrak{M}_{\text{0D}}$ for all 72 geometric models with Windkessel BCs
$\boldsymbol{\theta}_{\text{VMR}}$ from the Vascular Model Repository to
generate the ground truth data. We extracted the solution of pressures and
flow rates and their time derivatives from each model at all node points of
the 0D LPNs. We initialized all 0D parameters $\boldsymbol{\alpha}=\mathbf{0}$
to not introduce any bias to the optimized parameters
$\hat{\boldsymbol{\alpha}}$. Providing
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\boldsymbol{\alpha}_{\text{b}})$
as input, we performed LM optimization and obtained results for the optimized
branch parameters $\hat{\boldsymbol{\alpha}}_{\text{b}}$.

Due to the efficient C++ implementation and fast Python interface of
svZeroDSolver, running all 72 forward 0D simulations and optimizations took
only 47 s and 101 s on a single core, respectively. The optimization typically
converged in less than ten iterations, although three coronary models
stagnated at residual norms slightly above the tolerances defined in Appendix
B.2. This could be due to pressures and flow rates spanning different orders
of magnitude in the aorta and coronary arteries, which could be prevented by
normalization. We provide a detailed overview of all optimized parameters
$\hat{\boldsymbol{\alpha}}_{\text{b}}$ vs. ground truth geometric parameters
$\boldsymbol{\alpha}_{\text{b}}$ in Figure 12 in Appendix D, showing
$\mathcal{R},\mathcal{L},\mathcal{C},$ and $\mathcal{S}$ for every branch. The
mean of the coefficients of determination averaged over all 72 models and
their branches was
$R^{2}_{\mathcal{R}}=R^{2}_{\mathcal{L}},=R^{2}_{\mathcal{C}}=1.0$. For the
stenosis coefficient, this value was slightly reduced to
$R^{2}_{\mathcal{S}}=0.95$.

These results confirm that the four parameters of the BloodVessel element,
$\mathcal{R},\mathcal{L},\mathcal{C}$, and $\mathcal{S}$ are identifiable from
the solution vector $\mathbf{y}$ and its time derivative $\dot{\mathbf{y}}$.
The slightly lower value $R^{2}_{\mathcal{S}}$ of the stenosis coefficient
indicates that it is more difficult to identify. However, it should be noted
that many of the stenosis parameters $\mathcal{S}$ in the geometric 0D models
are negligibly small and have little influence on the result, which makes them
difficult to identify. The stenosis parameter is well identifiable in models
with pathological stenoses, like in aortic coarctations.

####  3.1.2 Optimizing zero-dimensional models to three-dimensional data

We demonstrated in 3.1.1 the excellent identifiability of 0D parameters from
ground truth data using the LM optimization. In this section, we quantify how
well optimized 0D models $\mathfrak{M}_{\text{0D}}$ approximate high-fidelity
3D models $\mathfrak{M}_{\text{3D}}$:

$\displaystyle\text{Compare\leavevmode\nobreak\
}\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{VMR}})\text{\leavevmode\nobreak\
vs.\leavevmode\nobreak\
}\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})\text{\leavevmode\nobreak\
with\leavevmode\nobreak\
}\hat{\boldsymbol{\alpha}}=\arg\min_{\boldsymbol{\alpha}}\|\mathbf{r}(\boldsymbol{\alpha},\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{VMR}}))\|^{2}.$
(27)

We extracted $\mathfrak{M}_{\text{3D}}$ velocity and pressure from the last
cardiac cycle at all time steps with Windkessel BC parameters
$\boldsymbol{\theta}_{\text{VMR}}$ given from the Vascular Model Repository.
Using a previously proposed technique to extract 3D results along the
centerline of the blood vessels, we mapped the cross-section averaged flow
rate and pressure from 3D to 0D LPN node points [3, 42]. We initialized the 0D
parameters with the values from the geometric model as
$\boldsymbol{\alpha}^{0}_{\text{b}}=\boldsymbol{\alpha}_{\text{b}}$ and
$\boldsymbol{\alpha}^{0}_{\text{j}}=\mathbf{0}$ in branches and junctions,
respectively. For model 0081_0001, we excluded the stenosis coefficient from
LM optimization as the 0D model with optimized stenosis coefficient did not
converge during forward evaluations. The optimization of all 72 models took 79
s in total on a single core.

Figure 13 in Appendix D provides a detailed correlation of all optimized 0D
parameters $\hat{\boldsymbol{\alpha}}_{\text{b}}$ vs. the geometric ones
$\boldsymbol{\alpha}_{\text{b}}$. The mean of the coefficients of
determination averaged over all 72 models were
$R^{2}_{\mathcal{R}}=0.68,R^{2}_{\mathcal{L}}=0.89,R^{2}_{\mathcal{C}}=0.23$,
and $R^{2}_{\mathcal{S}}=0.24$. This optimization also incorporates 0D
junction parameters $\hat{\boldsymbol{\alpha}}_{\text{j}}$. However, since
$\boldsymbol{\alpha}_{\text{j}}=\mathbf{0}$ in the geometric 0D models, we
cannot provide a correlation for these values. Since we aim to minimize the 0D
approximation errors by optimizing element parameters (27), we do not expect
the coefficients of determination to be close to one as in Section 3.1.1.
Rather, they indicate where the model generation process from a 3D geometry
can still be improved. Comparing the optimized parameters
$\hat{\boldsymbol{\alpha}}_{\text{b}}$ to the geometric ones
$\boldsymbol{\alpha}_{\text{b}}$ in Figure 13, there are some systematic
differences. As expected, due to the rigid wall nature of the 3D results, the
capacitances $\mathcal{C}$ are optimized as near zero. Furthermore, the
optimized stenosis parameters $\mathcal{S}$ are consistently higher,
indicating that they are systematically underestimated in the geometric
models. However, some pulmonary models with higher post-optimization errors
than others exhibit non-zero capacitances $\mathcal{C}$. Given that these
models have many more junctions than others, this could indicate limitations
in our BloodVesselJunction model in representing 3D fluid dynamics in a
junction.

Figure 3 depicts the change in the errors between geometric 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\boldsymbol{\alpha}_{\text{b}})$
and optimized 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})$
with respect to 3D models
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{VMR}})$. The optimization
improved the maximum approximation errors for pressures
$\epsilon_{P,\text{max}}$ and flow rates $\epsilon_{Q,\text{max}}$ in 71 out
of the 72 models, with error metrics defined as [3]

$\epsilon_{P,\text{max}}=\frac{n_{t}}{n_{\text{cap}}}\sum_{i=1}^{n_{\text{cap}}}\frac{\max_{t}\left|P_{t,i}^{\text{0D}}-P_{t,i}^{3\text{D}}\right|}{\sum_{t=1}^{n_{t}}P_{t,i}^{\text{3D}}},\quad\epsilon_{Q,\text{max}}=\frac{1}{n_{\text{out}}}\sum_{i=1}^{n_{\text{out}}}\frac{\max_{t}\left|Q_{t,i}^{\text{0D}}-Q_{t,i}^{3\text{D}}\right|}{\max_{t}Q_{t,i}^{\text{3D}}-\min_{t}Q_{t,i}^{\text{3D}}},\\\
$ (28)

with $n_{t}$ time steps $t$, number of outlets $n_{\text{out}}$, and number of
caps $n_{\text{cap}}=n_{\text{out}}+1$. Only the coronary model 0186_0002
performed slightly worse after optimization. For this model, the junction
between the aorta and the coronary arteries could not be calibrated properly.
Overall, the median error was reduced from 5.6% to 0.6% for pressure and 12.3%
to 1.4% for flow following optimization. The optimization achieved the largest
improvements for aortofemoral models. Given that our approach optimizes all 0D
parameters, it achieves much lower approximation errors than previous studies
where only resistances were adapted, e.g., to match key 3D outputs within
$10\%$ [13].

Optimizing the 0D models to 3D data lowered their maximum approximation error
for pressure and flow rate by nearly one order of magnitude, from $\sim 10\%$
to $\sim 1\%$. This considerable improvement in 0D accuracy is key to our
Bayesian Windkessel calibration method. The improvements can be mainly
attributed to the optimization of previously unknown BloodVesselJunction
parameters $\hat{\boldsymbol{\alpha}}_{\text{j}}$. Since the geometric model
does not capture the fluid behavior in the junctions, it consistently
underestimates their pressure drop. Informing the 0D model using a 3D
simulation result, it is easily possible to also derive suitable parameters
for the junctions and thereby improve the 0D model to model the full domain
accurately.

Figure 3: Changes in maximum pressure and flow error between geometric and
optimized 0D models. Blue is an improvement; red is a deterioration. Darker
values indicate larger changes. Statistical results are provided on the right.
See Figure 7 for an overview of the geometries.

####  3.1.3 Generalization of optimized zero-dimensional models

We demonstrated in Section 3.1.2 that the LM optimization yields 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})$
which reduce maximum pressure and flow errors with respect to 3D models
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{VMR}})$ by nearly one
order of magnitude compared to geometric models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\boldsymbol{\alpha}_{\text{b}})$
for fixed Windkessel parameters $\boldsymbol{\theta}_{\text{VMR}}$. In this
section, we demonstrate that 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{v}},\hat{\boldsymbol{\alpha}})$
optimized with Windkessel training parameters $\boldsymbol{\theta}_{\text{t}}$
generalize well to a wide range of Windkessel validation parameters
$\boldsymbol{\theta}_{\text{v}}$:

$\displaystyle\text{Compare\leavevmode\nobreak\
}\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{v}})\text{\leavevmode\nobreak\
vs.\leavevmode\nobreak\
}\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{v}},\hat{\boldsymbol{\alpha}})\text{\leavevmode\nobreak\
with\leavevmode\nobreak\
}\hat{\boldsymbol{\alpha}}=\arg\min_{\boldsymbol{\alpha}}\|\mathbf{r}(\boldsymbol{\alpha},\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{\text{t}}))\|^{2},\quad\boldsymbol{\theta}_{\text{v}}\neq\boldsymbol{\theta}_{\text{t}}.$
(29)

This generalization is crucial to our goal of using 0D models as a surrogate
in Bayesian Windkessel calibration. If optimized 0D models generalize well
over a wide range of Windkessel parameters $\boldsymbol{\theta}$, we can
confidently approximate the posterior distribution
$\mathfrak{p}_{\text{3D}}\approx\mathfrak{p}_{\text{0D}}$.

Since this generalization study requires extensive sampling of computationally
expensive 3D simulation results, we select a subset of three simulation
models: the aortic model of an 11-year-old human female with aortic
coarctation 0104_0001, the aortofemoral model of a 76-year-old human male with
an abdominal aortic aneurysm 0140_2001, and the pulmonary model of a 43-year-
old healthy human female 0080_0001. See Figure 7 for their geometries. These
models span three different anatomies and exhibit high pressure and flow
approximation errors of the geometric 0D models of around 10%, highlighted
with arrows in Figure 3. For each model, we reused a set
$\\{\boldsymbol{\theta}_{i}\\}$ of 50 variations of BC parameters from
Pegolotti et al. [45]. This dataset was generated by multiplying flow rates,
resistance, and capacitance values of the BCs parameters by a factor randomly
and independently sampled from a uniform distribution between 0.8 and 1.2. In
each comparison, we used identical BCs parameters $\boldsymbol{\theta}_{i}$ in
geometric 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{i},\boldsymbol{\alpha}_{\text{b}})$
and optimized ones
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{i},\hat{\boldsymbol{\alpha}})$
and compared them w.r.t. the corresponding 3D models
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta}_{i})$. To ensure comparability,
we applied the 3D result of the first time step in the last simulated cardiac
cycle as initial conditions to the 0D simulation. The 0D models were then
evaluated for one cardiac cycle to ensure that the cardiac cycles in 0D and 3D
were identical.

For each validation, we optimized the 0D model to one BC training parameter
set $\boldsymbol{\theta}_{\text{t}}$ and compared approximation errors for all
other 49 BC validation parameters $\boldsymbol{\theta}_{\text{v}}$:

$\displaystyle\boldsymbol{\theta}_{\text{t},i}=\boldsymbol{\theta}_{i},\quad\boldsymbol{\theta}_{\text{v},i}=\\{\boldsymbol{\theta}_{1},\dots,\boldsymbol{\theta}_{50}\\}\setminus\boldsymbol{\theta}_{i},\quad\forall
i\in\\{1,\dots,50\\}.$ (30)

We perform this validation fifty times to obtain a cross-validation. Figure 4
shows the cross-validation error at systole for pressure and flow. The
performance of the optimized model is further divided into performance on the
training set $\\{\boldsymbol{\theta}_{\text{t},i}\\}$ and the validation set
$\\{\boldsymbol{\theta}_{\text{v},i}\\}$. For reference, we also show errors
for the geometric 0D model, i.e., without optimization. There is a large
variance of errors for the geometric model evaluated on all
$|\\{\boldsymbol{\theta}_{i}\\}|=50$ BCs parameter sets. As observed in
Section 3.1.2, the error is drastically reduced on the training set of all
three optimized 0D models. Furthermore, the variance of the
$|\\{\boldsymbol{\theta}_{\text{t},i}\\}|=50$ training errors is much lower.
Most importantly, the validation error is equal to or slightly larger than the
training error. Despite encompassing
$|\\{\boldsymbol{\theta}_{\text{v},i}\\}|=50\times 49=2450$ simulations, the
variance of the validation error is only slightly larger than the one of the
50 training errors.

We summarize two observations: First, although the optimized pulmonary models
have been significantly improved, they perform worse than the aortic models.
This could be related to the increased size and, thus, the complexity of the
model but could also hint toward a limitation of the junction model used
herein, as junctions make up a large portion of the pulmonary model domain.
Second, we do not observe a notable difference in the low 0D model error from
Section 3.1.2 even for randomly sampled BC parameters. This confirms good
generalization of the LM optimization approach, which demonstrates that we can
confidently use the approximation
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})\approx\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$
for the Bayesian Windkessel BC calibration in Section 3.2.

Figure 4: Cross-validation of optimized 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})$
and comparison to geometric 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\boldsymbol{\alpha}_{\text{b}})$.
In each validation, the model was calibrated to the test parameters, and the
result was compared to the 49 validation parameters.

###  3.2 Bayesian Windkessel calibration incorporating optimized zero-
dimensional models

This section applies the BC calibration approach proposed in Section 2.1 to
the aortic model of an 11-year-old human female with aortic coarctation
0104_0001. We split this section into problem-setup (3.2.1), results (3.2.2),
and validation (3.2.3).

####  3.2.1 Problem setup

Our problem setup is typical for cardiovascular fluid dynamics simulations and
is often performed manually or according to well-known heuristics. Typical
flow measurements are obtained from 4D Flow Magnetic Resonance Imaging at the
inlet and outlets of the model [46]. Pressure measurements at the inlet are
obtained either from in vivo catheter measurements or estimated from cuff
measurements [47].

In this section, one BC parameter denoted as $\theta^{(i)}$ is calibrated for
each outlet (i.e., each Windkessel) of the model. The parameter $\theta^{(i)}$
is the logarithmic total resistance at the respective Windkessel outlet

$\theta^{(i)}=\log(R^{(i)}),\quad
R^{(i)}=R_{\text{p}}^{(i)}+R_{\text{d}}^{(i)}.$ (31)

The logarithmic reparameterization ensures that the resistances remain
positive, independent of the parameter choice, as negative resistances are
unphysiological and usually lead to failed forward model evaluations. This
yields the overall parameter vector

$\boldsymbol{\theta}=\left[\begin{matrix}\log\left(R^{(i)}\right),&\dots,&\log\left(R^{(n_{\text{out}})}\right)\end{matrix}\right]^{\top},$
(32)

with the number of outlets $n_{\text{out}}$. The individual parameters of each
Windkessel model, i.e., $R_{\text{p}}$, $R_{\text{d}}$, and $C$, are derived
from the total resistance based on two assumptions: The ratio between distal
to proximal resistance is assumed as constant and given

$\frac{R_{\text{d},i}}{R_{\text{p},i}}=\mathrm{const},$ (33)

and the time constant $\tau_{i}$ of each Windkessel BC is known and constant

$\tau_{i}=R_{\text{d},i}\cdot C_{i}=\mathrm{const}.$ (34)

The respective values for $\tau_{i}$ and
$\frac{R_{\text{d},i}}{R_{\text{p},i}}$ are unique for each Windkessel BC and
are calculated based on the existing BC parameters
$\boldsymbol{\theta}_{\text{VMR}}$ in the Vascular Model Repository. The
observations $\boldsymbol{y}_{\text{obs}}$ for this example are artificially
created based on one 3D simulation for the geometry with BC parameters
$\boldsymbol{\theta}_{\text{VMR}}$. In agreement with typical clinical
measurements, the following observations have been selected: the minimum
pressure at the inlet, the maximum pressure at the inlet, and the time-
averaged flow rates at each outlet

$\boldsymbol{y}_{\text{obs}}=\left[\begin{matrix}P_{\text{in,min}},&P_{\text{in,max}},&Q_{\text{mean,1}},&\dots,&Q_{\text{mean},n_{\text{outlets}}}\end{matrix}\right]^{\top}.$
(35)

The aortic anatomy 0104_0001 has five outlets connected to Windkessel BCs and
accordingly five calibration parameters $\theta^{(i)}$ and seven observations
$y_{\text{obs}}^{(i)}$. It is assumed that the observations are corrupted by
additive Gaussian noise as specified in Section 2.1.2, which is characterized
by a constant signal-to-noise ratio

$\text{SNR}_{i}=\frac{y_{\text{obs},i}^{2}}{\sigma_{i}^{2}}=\text{const},$
(36)

for all observations $y_{\text{obs},i}$.

####  3.2.2 Posterior distribution for Windkessel boundary conditions

The Bayesian Windkessel calibration is performed with three signal-to-noise
ratios of 100, 11.1, and 4 to investigate changes in the posterior
distribution for increasing noise in observations
$\boldsymbol{y}_{\text{obs}}$. SMC is performed with $k=10\,000$ particles, a
resampling threshold $\text{ESS}_{\text{min}}=5\,000$, and two rejuvenation
steps per iteration. All calculations were performed on Stanford’s high
performance computing cluster Sherlock using four 12-core Intel Xeon Gold 5118
CPUs. It should be noted, however, that only the 3D simulation took advantage
of all 48 available cores. We ran the SMC algorithm on 24 cores and the LM
calibration on a single core. The parameter estimation for the
$\text{SNR}=100$ case took a bit over 24 h, where 1.5h was spent on the first
SMC run (with 210 000 0D evaluations), 21.5 h on the 3D simulation and 1.5 h
on the second SMC run. The LM calibration took only 20 s.

Figure 5 presents the posterior densities for all parameters $\theta^{(i)}$
and signal-to-noise ratios SNR. All parameter distributions are concentrated
around single values for the high signal-to-noise ratio $\text{SNR}=100$. The
locations of the peaks match well with the parameter values used to create the
artificial observations (marked by the crossing gray lines). The posterior
variance is low, which matches the expectation given that the noise in
observations is low. The posterior variance is generally larger for the lower
single-to-noise ratio $\text{SNR}=11.1$, and some distributions exhibit multi-
modality. The posterior distributions for the BC parameter at the descending
aorta $\theta^{(2)}$ are narrower compared to the other parameters, suggesting
a higher sensitivity of the model output to $\theta^{(2)}$. For the largest
noise case $\text{SNR}=4$, the posterior is characterized by the largest
variance.

The results show that the proposed Bayesian calibration approach successfully
captured the multi-dimensional posterior distribution of the Windkessel
parameters based on the optimized 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta}_{\text{VMR}},\hat{\boldsymbol{\alpha}})$
and delivered robust results even for completely uninformed and broad priors.
It has been demonstrated that the Bayesian calibration approach is
advantageous in settings where observations are corrupted with noise. While
the posterior collapses almost to a Dirac distribution in the case of a high
signal-to-noise ratio in the observations, the cases with high measurement
noise highlight how the uncertainty propagates to the posterior.

Calculating the full posterior distribution provides many advantages compared
to deterministic optimization strategies where only a point estimate is
obtained. First, the probabilistic calibration problem is better posed and
less prone to get stuck in local optima. Second the posterior distribution
provides a quantification of confidence in the calibration result. This is
instrumental for future clinical application of these models as it allows one
to assess the expected reliability of a simulation model calibrated based on
very sparse data. Knowing the full posterior distribution allows for more
nuanced and better-informed parameter choices. Third, the posterior
distributions can be used to quantify uncertainty in model predictions using
forward propagation, e.g., wall shear stress in thoracic aortic aneurysms [48]
or pressure drop over an aortic coarctation [49]. In the latter, the
uncertainty in the model prediction can be taken into account while making the
clinical decision of whether to intervene or not [4]. We enabled the sampling
intensive SMC algorithm by developing a fast svZeroDSolver that performed the
hundreds of thousand 0D evaluations in a bit over an hour.

Figure 5: 2D marginals of the 5D posterior for all parameters $\theta^{(i)}$
and signal-to-noise ratios SNR based on model 0104_0001. The diagonal plots
show the 1D kernel density estimate. The upper right plots the particles
colored by their weight, and the lower left plots the 2D density contours. The
parameter values corresponding to the ground truths are marked with crossing
gray lines.

####  3.2.3 Comparison of 0D posterior with 3D posterior

The framework proposed in Section 2.1 relies on the assumption that 0D models
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$ with
optimized 0D parameters $\hat{\boldsymbol{\alpha}}$ approximate the 3D model
$\mathfrak{M}_{\text{3D}}(\boldsymbol{\theta})$ accurately. This was
demonstrated in a deterministic sense in Section 3.1 for a wide range of
Windkessel BC parameters $\boldsymbol{\theta}$. In this section, we validate
that a posterior distribution based on an optimized 0D model
$\mathfrak{p}_{\text{0D}}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta},\hat{\boldsymbol{\alpha}})$
agrees with the posterior distributions based on the full-dimensional 3D model
$\mathfrak{p}_{\text{3D}}(\boldsymbol{y}_{\text{obs}}|\boldsymbol{\theta})$
given identical uncertain observations $\boldsymbol{y}_{\text{obs}}$.

To ensure computational feasibility, we reduce the problem setup from Sections
3.2.1 and 3.2.2 to two dimensions. The parameters $\theta^{(1)}$,
$\theta^{(2)}$, $\theta^{(4)}$ and $\theta^{(5)}$ are now coupled and denoted
as $\theta^{(1,2,4,5)}$ in the following and $\theta^{(3)}$ is calibrated
individually. We evaluate the unnormalized posterior point-wise on a uniform
10x10 grid in the interval $[2,8]$ for $\theta^{(1,2,4,5)}$ and $[3,6]$ for
$\theta^{(3)}$ on a geometric 0D model
$\mathfrak{M}_{\text{0D}}(\boldsymbol{\alpha}_{\text{b}})$, an optimized 0D
model $\mathfrak{M}_{\text{0D}}(\hat{\boldsymbol{\alpha}})$, and a full-
dimensional 3D model $\mathfrak{M}_{\text{3D}}$. Afterward, we normalize the
unnormalized posterior to yield an estimate of the true posterior
distribution. A signal-to-noise ratio of $\text{SNR}=11.1\%$ is used.
Performing all 100 3D simulations took nearly 800 hours on two 64-core AMD
EPYC 7742 CPUs, while the 100 0D simulations were completed in less than a
minute on a single core.

Figure 6: Comparison of the posterior distribution evaluated on the geometric
0D model, the calibrated 0D model, and the 3D model. The uniform prior is
added in orange to the marginal distributions.

The resulting posterior distributions for the three models are presented in
Figure 6. The geometric and optimized 0D models exhibit more posterior
variance than the 3D model. While the geometric 0D model posterior shows no
overlap with the 3D posterior, the optimized 0D model has a much improved
agreement with the latter. Despite the slight elongation of the posterior for
the 0D models, 3D and optimized 0D models attribute most of the probability
mass in a similar region close to the coordinate (6, 5). As seen in Section
3.1.2, model 0104_0001 was intentionally selected as one of the worst-
performing 0D models. However, it still achieved a good approximation of the
3D posterior, especially compared to the geometric 0D model, with a tendency
to overestimate the uncertainty in the calibration result.

The optimized 0D posterior deviates slightly from the high-fidelity 3D
posterior but clearly encapsulates the latter. In contrast, the geometric 0D
model completely fails to provide a viable surrogate for the 3D posterior.
While the geometric 0D model attributes the majority of its probability
density to wrong parameter choices and, thus, provides confidence in incorrect
results, the optimized 0D model merely overpredicts the posterior variance.
Informing the 0D model by a single 3D simulation result was sufficient to
increase the agreement between 0D model and 3D model to the extent that the 0D
model can confidently be used for parameter calibration without requiring
further 3D evaluations. These results confirm that the optimized 0D models are
not only superior to the geometric 0D models for approximating the 3D
posterior, but are essential. Moreover, it has been confirmed that the
approach proposed in this work allows for accurate and efficient probabilistic
parameter calibration for Windkessel BCs.

##  4 Limitations and outlook

We demonstrated that a LM optimization approach can derive highly accurate 0D
models based on a 3D simulation of a single heartbeat. The optimized 0D models
showed great agreement with the physiological behavior of 3D models even
during substantial changes in the BCs. The optimization also proved effective
for automatically determining parameters for 0D junctions, which is a known
challenge for cardiovascular reduced-order models. The results provide
valuable insights regarding the potential to improve the automatic creation of
LPN models. This could be used to develop more accurate 0D element parameters
that are informed by both geometry and hemodynamics or as training data for
machine-learning-based approaches to 0D model generation [38].

Despite the significant reduction in 0D-3D approximation errors, there are
still some deviations between the optimized 0D models and the 3D models,
especially for pulmonary models. While it is not possible to achieve perfect
agreement between 3D and 0D models, this indicates that improvements to the 0D
models themselves may need to be investigated to include more nonlinear terms
that better capture 3D flow physics. A structured analysis of the
physiological effects currently not accounted for by the 0D elements could
result in better agreement. Additionally, a probabilistic optimization
approach could provide more information about the identifiability of the
individual 0D parameters based on the 3D training data. The 0D model
optimization has been performed using one cardiac cycle of a periodically
converged simulation. Future studies should investigate whether the amount of
3D data necessary to sufficiently inform the 0D model optimization can be
further reduced, i.e., a few 3D time steps may suffice. This would reduce the
computational expense of the algorithm even further.

We demonstrated that optimized 0D models can be deployed for accurate and
efficient Bayesian Windkessel calibration using SMC. The Bayesian approach not
only provides more robustness for the calibration, but the resulting posterior
distributions also provide a powerful tool for making more informed parameter
decisions in patient-specific cardiovascular modeling. While deterministic
calibration results can only be interpreted in one way (as the final parameter
choice), the posterior distributions can be processed to retrieve the optimal
parameter choice in view of multiple aspects (e.g. risk-assessment,
feasibility analysis). The posterior also provides additional information
about the uncertainty in the parameters given potentially noisy clinical
measurements. This is especially relevant in a clinical setting where data is
usually sparse and noisy. The results can also be used to decide whether a
parameter choice can be made confidently based on the given clinical
measurement or if more data is required.

We evaluated the proposed calibration framework on artificially created
measurement data. While the problem setup aimed to match a common calibration
scenario, further work could focus on using subject data. Applying the
proposed calibration framework to different setups can validate its robustness
on realistic clinical data. Moreover, we limited our investigations to rigid-
wall 3D models of the vasculature. While this is a valid assumption for the
stiffer vasculature or lower pressures, it is not generally applicable.
Despite this limitation, we previously demonstrated that our forward [3] and
inverse 0D models can incorporate vessel wall elasticity.

Furthermore, the posterior distribution derived based on the optimized 0D
models has been compared against the high-fidelity 3D posterior and proved to
be a valid surrogate for the latter. The poor performance of the geometric 0D
model highlighted that the 0D-3D optimization is essential for performing
Bayesian calibration on 0D models. Further work could generalize the algorithm
presented herein using an Expectation-Maximization approach that treats the 0D
model parameters $\boldsymbol{\alpha}$ as latent variables that are
consistently optimized as part of the inverse problem.

## Acknowledgements

This work was supported by grants K99HL161313, R01HL167516, and R01HL141712
from the National Heart, Lung, and Blood Institute and the Stanford Maternal
and Child Health Research Institute. DES acknowledges partial support from NSF
CAREER award #1942662 and NSF CDS&E award #2104831. WAW and JN acknowledge
support from the European Research Council (grant agreement No.
101021526-BREATHE). The authors gratefully acknowledge the Stanford Research
Computing Center and the San Diego Super Computing Center (SDSC) Expanse
cluster for providing the computational resources necessary for the numerical
simulations presented in this work.

##  Appendix A Overview of cardiovascular models

Figure 7 shows an overview of all 72 models used in this work, identical to
the ones in Pfaller et al. [3]. They are openly available from the Vascular
Model Repository [44].

Figure 7: Overview of models used in this work, sorted by anatomy. This set is
identical to the one in Pfaller et al. [3].

##  Appendix B Computational details of zero-dimensional models

This section contains the zero-dimensional (0D) time integration scheme (B.1)
and the system matrices for the BloodVessel (B.3), BloodVesselJunction (B.4),
and Windkessel (B.5) lumped-parameter network elements as implemented in
svZeroDSolver [41].

###  B.1 Forward time integration

We use a generalized-$\alpha$ time integration [50, 51] and solve for the
solution vector $\mathbf{y}$ in each time step $n$ using the Newton-Raphson
method,

$\displaystyle\text{Predict:\leavevmode\nobreak\
}\quad\dot{\mathbf{y}}_{n+1}^{0}$
$\displaystyle=\frac{\gamma-1}{\gamma}\dot{\mathbf{y}}_{n},$ (37)
$\displaystyle\mathbf{y}_{n+1}^{0}$ $\displaystyle=\mathbf{y}_{n},$ (38)
$\displaystyle\text{Initiate:\leavevmode\nobreak\
}\quad\dot{\mathbf{y}}_{n+\alpha_{m}}^{i}$
$\displaystyle=\dot{\mathbf{y}}_{n}+\alpha_{m}\left(\dot{\mathbf{y}}_{n+1}^{i}-\dot{\mathbf{y}}_{n}\right),$
(39) $\displaystyle\mathbf{y}_{n+\alpha_{f}}^{i}$
$\displaystyle=\mathbf{y}_{n}+\alpha_{f}\left(\mathbf{y}_{n+1}^{i}-\mathbf{y}_{t}\right),$
(40) $\displaystyle\text{Solve:\leavevmode\nobreak\
}\quad\mathbf{K}(\mathbf{y}_{n+\alpha_{f}}^{i},\dot{\mathbf{y}}_{n+\alpha_{m}}^{i})\cdot\Delta\dot{\mathbf{y}}_{n+1}^{i}$
$\displaystyle=-\mathbf{r}(\mathbf{y}_{n+\alpha_{f}}^{i},\dot{\mathbf{y}}_{n+\alpha_{m}}^{i}),$
(41) $\displaystyle\text{Update:\leavevmode\nobreak\
}\quad\dot{\mathbf{y}}_{n+1}^{i+1}$
$\displaystyle=\dot{\mathbf{y}}_{n+1}^{i}+\Delta\dot{\mathbf{y}}_{n+1}^{i},$
(42) $\displaystyle\mathbf{y}_{n+1}^{i+1}$
$\displaystyle=\mathbf{y}_{n+1}^{i}+\Delta\dot{\mathbf{y}}_{n+1}^{i}\,\gamma\Delta
t_{n},$ (43)

with time step $\Delta t_{n}$, residual $\mathbf{r}$, and Jacobian
$\mathbf{K}$. We derive the generalized-$\alpha$ parameters $\alpha_{f}$,
$\alpha_{m}$, and $\gamma$ from the spectral radius $\rho_{\infty}$

$\displaystyle\alpha_{m}=\frac{1-\rho_{\infty}}{2(1+\rho_{\infty})},\quad\alpha_{f}=\frac{1}{1+\rho_{\infty}},\quad\gamma$
$\displaystyle=\frac{1}{2}+\alpha_{m}-\alpha_{f}.$ (44)

###  B.2 Inverse parameter optimization

We use the Levenberg-Marquardt (LM) algorithm [52, 53] to solve the least
squares problem (24) iteratively for the optimal parameter set
$\hat{\boldsymbol{\alpha}}$,

$\displaystyle\text{update\leavevmode\nobreak\ }\boldsymbol{\alpha}^{i+1}$
$\displaystyle=\boldsymbol{\alpha}^{i}+\Delta\boldsymbol{\alpha}^{i+1},$ (45)
$\displaystyle\text{solve\leavevmode\nobreak\
}\left[\mathbf{J}^{\mathrm{T}}\mathbf{J}+\lambda\operatorname{diag}\left(\mathbf{J}^{\mathrm{T}}\mathbf{J}\right)\right]^{i}\cdot\Delta\boldsymbol{\alpha}^{i+1}$
$\displaystyle=-\left[\mathbf{J}^{\mathrm{T}}\mathbf{r}\right]^{i},\quad\lambda^{i}=\lambda^{i-1}\cdot\left\|\left[\mathbf{J}^{\mathrm{T}}\mathbf{r}\right]^{i}\right\|_{2}/\left\|\left[\mathbf{J}^{\mathrm{T}}\mathbf{r}\right]^{i-1}\right\|_{2}.$
(46)

The LM algorithm approximates the Hessian from the Jacobian matrix
$\mathbf{J}$ as $\Delta S\approx\mathbf{J}^{\mathrm{T}}\mathbf{J}$. The second
term introduces damping with factor $\lambda$, which we update based on the
change in gradient $\nabla S=\mathbf{J}^{\mathrm{T}}\mathbf{r}$. Thus, the
damping parameter approaches $\lambda\to 0$ as the parameter set
$\boldsymbol{\alpha}\to\hat{\boldsymbol{\alpha}}$ approaches the optimal
solution. For $\lambda\to\infty$, we obtain the steepest descent method, and
for $\lambda=0$, we obtain the Gauss-Newton method. We terminate the
optimization (46) when both $L^{2}$ norms,

$\displaystyle\left\|\left[\mathbf{J}^{\mathrm{T}}\mathbf{r}\right]^{\mathrm{i}}\right\|_{2}<\operatorname{tol}_{\text{grad}}^{\boldsymbol{\alpha}}\text{
and
}\left\|\Delta\boldsymbol{\alpha}^{\mathrm{i}+1}\right\|_{2}<\operatorname{tol}_{\text{inc}}^{\boldsymbol{\alpha}},$
(47)

are below the tolerances for gradient and increment,
$\operatorname{tol}_{\text{grad}}^{\boldsymbol{\alpha}}$ and
$\operatorname{tol}_{\text{inc}}^{\boldsymbol{\alpha}}$, respectively. In this
work, we chose the LM optimization parameters of initial damping factor
$\lambda^{0}=1$, gradient tolerance
$\operatorname{tol}_{\text{grad}}^{\boldsymbol{\alpha}}=10^{-5}$, increment
tolerance $\operatorname{tol}_{\text{inc}}^{\boldsymbol{\alpha}}=10^{-10}$,
and a maximum number of iterations of $100$.

###  B.3 BloodVessel element

We model the blood vessel branches with the BloodVessel element, whose circuit
is visualized in Figure 8.

$Q_{\text{in}}$$P_{\text{in}}$$\mathcal{R}$$\mathcal{S}$$\mathcal{L}$$P_{\text{out}}$$\mathcal{C}$$Q_{\text{out}}$
Figure 8: The BloodVessel element with resistance $\mathcal{R}$, stenosis
$\mathcal{S}$, inertance $\mathcal{L}$, and capacitance $\mathcal{C}$.

The governing equations are given by

$\displaystyle-(\mathcal{R}+\mathcal{S}|Q_{\text{in}}|)\,Q_{\text{in}}-\mathcal{L}\,\dot{Q}_{\text{out}}-\Delta
P$ $\displaystyle=0,$ (48)
$\displaystyle\mathcal{C}(\mathcal{R}+2\mathcal{S}|Q_{\text{in}}|)\,\dot{Q}_{\text{in}}-\mathcal{C}\dot{P}_{\text{in}}-\Delta
Q$ $\displaystyle=0.$ (49)

The remainder of this section contains the forward (B.3.1) and inverse (B.3.2)
problem of the BloodVessel element.

####  B.3.1 Forward problem

The vector of element unknowns is

$\displaystyle\mathbf{y}^{(e)}$
$\displaystyle=\left[\begin{array}[]{llll}P_{\text{in}}^{(e)}&Q_{\text{in}}^{(e)}&P_{\text{out}}^{(e)}&Q_{\text{out}}^{(e)}\end{array}\right]^{\mathrm{T}}$
(51)

The BloodVessel element matrices for the forward problem (18) are

$\displaystyle\mathbf{E}^{(e)}$
$\displaystyle=\left[\begin{array}[]{cccc}0&0&0&-\mathcal{L}\\\
-\mathcal{C}&\mathcal{C}\mathcal{R}&0&0\end{array}\right]^{(e)},$ (54)
$\displaystyle\mathbf{F}^{(e)}$
$\displaystyle=\left[\begin{array}[]{cccc}1&-\mathcal{R}&-1&0\\\
0&1&0&-1\end{array}\right]^{(e)},$ (57) $\displaystyle\mathbf{c}^{(e)}$
$\displaystyle=\mathcal{S}\,|Q_{\text{in}}|\left[\begin{array}[]{c}-Q_{\text{in}}\\\
2\mathcal{C}\dot{Q}_{\text{in}}\end{array}\right]^{(e)}.$ (60)

Only the stenosis term $\mathcal{S}$ contributes to the nonlinear term
$\mathbf{c}$. Its linearizations are

$\displaystyle\left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{(e)}$
$\displaystyle=\mathcal{S}\,|Q_{\text{in}}|\left[\begin{array}[]{cccc}0&0&0&0\\\
0&2\mathcal{C}&0&0\end{array}\right]^{(e)},$ (63)
$\displaystyle\left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{(e)}$
$\displaystyle=\mathcal{S}\,\text{sgn}(Q_{\text{in}})\left[\begin{array}[]{cccc}0&-2Q_{\text{in}}&0&0\\\
0&2\mathcal{C}\dot{Q}_{\text{in}}&0&0\end{array}\right]^{(e)}.$ (66)

####  B.3.2 Inverse problem

The vector of element unknowns is

$\displaystyle\boldsymbol{\alpha}^{(e)}$
$\displaystyle=\left[\begin{array}[]{llll}\mathcal{R}^{(e)}&\mathcal{C}^{(e)}&\mathcal{L}^{(e)}&\mathcal{S}^{(e)}\end{array}\right]^{\mathrm{T}},$
(68)

containing flow and pressure for the inlet and $n$ outlets. The BloodVessel
element Jacobian matrix for the inverse problem (23) is

$\displaystyle\mathbf{J}^{(e)}=\left[\begin{array}[]{cccc}-Q_{\text{in}}&0&-\dot{Q}_{\text{out}}&-|Q_{\text{in}}|Q_{\text{in}}\\\
\mathcal{C}\,\dot{Q}_{\text{in}}&-\dot{P}_{\text{in}}+(\mathcal{R}+2\mathcal{S}\,|Q_{\text{in}}|)\,\dot{Q}_{\text{in}}&0&2\mathcal{C}|Q_{\text{in}}|\dot{Q}_{\text{in}}\end{array}\right]^{(e)}.$
(71)

###  B.4 BloodVesselJunction element

We model the junctions between blood vessel branches with the
BloodVesselJunction element. It is visualized in Figure 9 (top) for an inlet
connecting to $1,\dots,n$ outlets. An individual junction model connecting the
inlet to outlet $i$ is visualized in Figure 9 (bottom).

$Q_{\text{in}}$$P_{\text{in}}$$BV_{1}$$P_{\text{out},1}$$BV_{n}$$P_{\text{out},n}$$Q_{\text{out},1}$$Q_{\text{out},n}$…

$Q_{\text{in}}$$P_{\text{in}}$$\mathcal{R}_{i}$$\mathcal{S}_{i}$$\mathcal{L}_{i}$$P_{\text{out},i}$$Q_{\text{out},i}$

Figure 9: The BloodVesselJunction element (top), connecting an inlet to $n$
outlets. An individual inlet-outlet connection $BV_{i}$ is shown at the
bottom, with resistance $\mathcal{R}_{i}$, stenosis $\mathcal{S}_{i}$, and
inertance $\mathcal{L}_{i}$.

The governing equations for the BloodVesselJunction are

$\displaystyle Q_{\text{in}}-\sum_{i}^{n}Q_{\text{out},i}$ $\displaystyle=0,$
(72) $\displaystyle\Delta
P_{i}+(\mathcal{R}+\mathcal{S}|Q_{\text{out},i}|)\,Q_{\text{out},i}+\mathcal{L}\,\dot{Q}_{\text{out},i}$
$\displaystyle=0,\quad\forall\,i=1,\dots,n.$ (73)

The remainder of this section contains the forward (B.4.1) and inverse (B.4.2)
problem of the BloodVesselJunction element.

####  B.4.1 Forward problem

The vector of element unknowns is

$\displaystyle\mathbf{y}^{(e)}$
$\displaystyle=\left[\begin{array}[]{lllllll}P_{\text{in}}^{(e)}&Q_{\text{in}}^{(e)}&P_{\text{out},1}^{(e)}&Q_{\text{out},1}^{(e)}&P_{\text{out},2}^{(e)}&Q_{\text{out},2}^{(e)}&\dots\end{array}\right]^{\mathrm{T}}.$
(75)

The BloodVesselJunction element matrices for the forward problem (18) are

$\displaystyle\mathbf{E}^{(e)}$
$\displaystyle=\left[\begin{array}[]{ccccccc}0&0&0&0&0&0&\dots\\\
0&0&0&-\mathcal{L}_{1}&0&0&\\\ 0&0&0&0&0&-\mathcal{L}_{2}&\\\
\vdots&&&&&&\ddots\end{array}\right]^{(e)},$ (80)
$\displaystyle\mathbf{F}^{(e)}$
$\displaystyle=\left[\begin{array}[]{ccccccc}0&1&0&-1&0&-1&\dots\\\
1&0&-1&-\mathcal{R}_{1}&0&0&\\\ 1&0&0&0&-1&-\mathcal{R}_{2}&\\\
\vdots&&&&&&\ddots\end{array}\right]^{(e)},$ (85)
$\displaystyle\mathbf{c}^{(e)}$ $\displaystyle=\left[\begin{array}[]{c}0\\\
-\mathcal{S}_{1}|Q_{\text{in},1}|\,Q_{\text{in},1}\\\
-\mathcal{S}_{2}|Q_{\text{in},2}|\,Q_{\text{in},2}\\\
\vdots\end{array}\right]^{(e)},$ (90)

with the first rows in each component containing the conservation of mass
equation (72). The linearization of the nonlinear term $\mathbf{c}$ is

$\displaystyle\left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{(e)}$
$\displaystyle=\mathbf{0},$ (91)
$\displaystyle\left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{(e)}$
$\displaystyle=2\left[\begin{array}[]{cccccccc}0&0&0&0&0&0&\dots\\\
0&0&0&-\mathcal{S}_{1}|Q_{\text{in},1}|&0&0&\\\
0&0&0&0&0&-\mathcal{S}_{1}|Q_{\text{in},2}|&\\\ \vdots&&&&&&\ddots\\\
\end{array}\right]^{(e)}.$ (96)

####  B.4.2 Inverse problem

The vector of element unknowns is

$\displaystyle\boldsymbol{\alpha}^{(e)}$
$\displaystyle=\left[\begin{array}[]{lllllllll}\mathcal{R}_{1}^{(e)}&\mathcal{R}_{2}^{(e)}&\dots&\mathcal{L}_{1}^{(e)}&\mathcal{L}_{2}^{(e)}&\dots&\mathcal{S}_{1}^{(e)}&\mathcal{S}_{2}^{(e)}&\dots\end{array}\right]^{\mathrm{T}},$
(98)

consisting of $3n$ unknowns for a junction with $n$ outlets. The
BloodVesselJunction element Jacobian matrix for the inverse problem (23) is

$\displaystyle\mathbf{J}^{(e)}=\left[\begin{array}[]{ccccccccc}0&0&\dots&0&0&\dots&0&0&\dots\\\
-Q_{\text{out},1}&0&&-\dot{Q}_{\text{out},1}&0&&|Q_{\text{in},1}|\,Q_{\text{in},1}&0&\\\
0&-Q_{\text{out},2}&&0&-\dot{Q}_{\text{out},2}&&0&|Q_{\text{in},2}|\,Q_{\text{in},2}&\\\
\vdots&&\ddots&\vdots&&\ddots&\vdots&&\ddots\end{array}\right]^{(e)}.$ (103)

###  B.5 Windkessel element

We use a three-element Windkessel visualized in Figure 10.

$Q_{\text{in}}$$P_{\text{in}}$$\mathcal{R}_{\text{p}}$$\mathcal{R}_{\text{d}}$$P_{\text{ref}}$$\mathcal{C}$
Figure 10: Three-element Windkessel BC

It is governed by the equation

$\displaystyle(\mathcal{R}_{\text{p}}+\mathcal{R}_{\text{d}})\,Q_{\text{in}}-P_{\text{in}}+P_{\text{ref}}-\mathcal{R}_{\text{d}}\,\mathcal{C}\,\dot{P}_{\text{in}}+\mathcal{R}_{\text{p}}\mathcal{R}_{\text{d}}\,\dot{Q}_{\text{in}}$
$\displaystyle=0.$ (104)

consisting of a capacitance $\mathcal{C}$, and proximal and distal resistances
$\mathcal{R}_{\text{p}}$ and $\mathcal{R}_{\text{d}}$, respectively [33]. In
our examples, the reference pressure is $P_{\text{ref}}=0$ at all outlets. The
vector of element unknowns for the Windkessel BC is

$\displaystyle\mathbf{y}^{(e)}$
$\displaystyle=\left[\begin{array}[]{ll}P_{\text{in}}^{(e)}&Q_{\text{in}}^{(e)}\end{array}\right]^{\mathrm{T}}.$
(106)

The forward problem (18) is given by

$\displaystyle\mathbf{E}^{(e)}$
$\displaystyle=\left[\begin{array}[]{cc}-\mathcal{R}_{\text{d}}\,\mathcal{C}&\mathcal{R}_{\text{p}}\,\mathcal{R}_{\text{d}}\end{array}\right]^{(e)},$
(108) $\displaystyle\mathbf{F}^{(e)}$
$\displaystyle=\left[\begin{array}[]{cccc}-1&\mathcal{R}_{\text{p}}+\mathcal{R}_{\text{d}}\end{array}\right]^{(e)},$
(110) $\displaystyle\mathbf{c}^{(e)}$
$\displaystyle=\left[\begin{array}[]{cccc}P_{\text{ref}}\end{array}\right]^{(e)}.$
(112)

The derivatives of the nonlinear term $\mathbf{c}$ vanish:

$\displaystyle\left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{(e)}=\left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{(e)}=\mathbf{0}.$
(113)

##  Appendix C Sequential Monte Carlo algorithm

Initialize:
$\boldsymbol{\theta}_{n}^{(i)}\sim\mathfrak{p}(\boldsymbol{\theta}),\,w_{n}^{(i)}=1,\,n=0,\,\gamma_{s}=0$$\pi_{n}=\left\\{\boldsymbol{\theta}_{n}^{(i)},\,w_{n}^{(i)}\right\\}$Evaluate:
$\pi_{n+1}^{\prime}=\mathfrak{p}\left(\boldsymbol{y}_{\text{obs}}|\mathfrak{M}(\boldsymbol{\theta}_{n}^{(1,\dots,k)})\right)$Tempering:
$\pi_{n+1}=(\pi_{n+1}^{\prime})^{\zeta_{n}},\,\gamma_{s}\mathrel{+}=\zeta_{n}$Reweight:
$w_{n+1}^{(i)}=w_{n}^{(i)}\frac{\pi_{n+1}^{(i)}}{\pi_{n}^{(i)}}$$\mathrm{ESS}<\mathrm{ESS_{\text{min}}}$Resample:
$\boldsymbol{\theta}_{n}^{(i)}=\mathcal{M}(W_{n}^{(1,\dots,k)})$Rejuvenate:
$\boldsymbol{\theta}_{n+1}^{(i)}=K\left(\boldsymbol{\theta}_{n}^{(i)}\right)$$\gamma_{s}\geq
1$Posterior:
$\mathfrak{p}(\boldsymbol{\theta}|y_{\text{obs}})\approx\sum_{i=1}^{k}W_{n+1}^{(i)}\,\delta_{\boldsymbol{\theta}_{n+1}^{(i)}}(\boldsymbol{\theta})$$n=n+1$YesNoYesNo
Figure 11: Sequential Monte Carlo algorithm.

Figure 11 gives an overview of the iterative Sequential Monte Carlo (SMC)
algorithm. Distributions in SMC are approximated by weighted particles
according to

$\pi_{n}(\boldsymbol{\theta})\approx\sum_{i=1}^{k}w_{n}^{(i)}\,\delta_{\boldsymbol{\theta}_{n}^{(i)}}(\boldsymbol{\theta}),$
(114)

wherein $k$ denotes the number of particles and $n$ is the iteration number.
Each particle $\delta_{\boldsymbol{\theta}_{n}^{(i)}}$ refers to a Dirac delta
distribution at particle location $\boldsymbol{\theta}_{i}$. Distributions
$\pi_{n}$ are not necessarily normalized and, as such, should not be regarded
as probability densities. A probability density is approximated with

$\mathfrak{p}_{n}(\boldsymbol{\theta})\approx\sum_{i=1}^{k}W_{n}^{(i)}\,\delta_{\boldsymbol{\theta}_{n}^{(i)}}(\boldsymbol{\theta})$
(115)

using the normalized weights

$W_{n}^{(i)}=\frac{w_{n}^{(i)}}{\sum_{j=1}^{k}w_{n}^{(j)}}.$ (116)

The algorithm is initialized by drawing $k$ particles from the prior
distribution with a uniform weight of one

$\boldsymbol{\theta}_{0}^{(i)}\sim\mathfrak{p}(\boldsymbol{\theta}),\,\,w_{0}^{(i)}=1,\quad\text{with}\quad
i\in\mathbb{N}_{>0},i\leq k.$ (117)

It is said that the particles target distribution $\pi_{0}$ (i.e., the prior).
By evaluating the likelihood of each particle, the successor distribution can
be derived

$\pi_{n+1}=\mathfrak{p}\left(y_{\mathrm{obs}}|\mathfrak{M}(\boldsymbol{\theta}_{n}^{(1,\dots,k)})\right)^{\zeta_{n+1}}.$
(118)

The particles are reweighted using importance sampling to target the new
distribution $\pi_{n+1}$

$w_{n+1}^{(i)}=w_{n}^{(i)}\cdot\frac{\pi_{n+1}^{(i)}}{\pi_{n}^{(i)}}.$ (119)

The exponent $\zeta_{n}$ refers to the tempering step size and is part of the
adaptive tempering reweighting strategy. It is selected based on the effective
sample size $\mathrm{ESS}$ and the resampling threshold
$\mathrm{ESS}_{\text{min}}$ according to

$\mathrm{ESS}(w_{n+1}^{(1,\dots,k)})\approx\zeta\mathrm{ESS}_{\text{min}}$
(120)

with

$\mathrm{ESS}=\frac{1}{\sum_{i=1}^{k}\left(w^{(i)}\right)^{2}}.$ (121)

The effective sample size measures how well the weights are distributed
between the particles and ensures an appropriate step size for each SMC step.
With every iteration, the weights of some particles increase further, and
others approach zero. Keeping many particles with (almost) no weights is
inefficient, so the particles are resampled. Particles with high weights are
duplicated, and particles with low weights are dropped. Herein, the systematic
resampling strategy from Chopin et al. [31] is adopted. The particles are
resampled based on a multinomial distribution
$\mathcal{M}(W_{n}^{(1,\dots,k)})$ when the effective sampling size
$\mathrm{ESS}$ is below the resampling threshold $\mathrm{ESS}_{\text{min}}$.
The probability of a particle ${\boldsymbol{\theta}^{(i)}}^{\prime}$ being
resampled from particle $\boldsymbol{\theta}^{(i)}$ is determined by its
relative weight

$\mathfrak{p}({\boldsymbol{\theta}^{(i)}}^{\prime})=\frac{w^{(i)}}{\sum_{j=1}^{k}w^{(j)}}=W^{(i)}.$
(122)

The resampling does not change the distribution targeted by the particles. To
avoid a high correlation between particles after resampling (as some particles
are duplicates of the same parameter vector), a rejuvenation step is added to
move all particles according to

$\boldsymbol{\theta}_{n+1}^{(i)}=K(\boldsymbol{\theta}_{n}^{(i)}),$ (123)

using transition kernel $K$ with stationary distribution $\pi_{n+1}$ according
to the proposal distribution in an MCMC step. The rejuvenation can be applied
multiple times in a row. The process is terminated when the tempering exponent
$\gamma_{s}$ reaches ones, defined as

$\gamma_{s}=\sum_{j=1}^{n+1}\zeta_{j}.$ (124)

##  Appendix D Detailed zero-dimensional optimization results

The results of the optimization of 0D parameters from 0D ground truth and
three-dimensional (3D) data are shown in Figures 12 and 13, respectively.

Figure 12: Correlation of optimized 0D elements from 0D ground truth data with
0D ground truth for all $N=72$ models (aorta, aortafemoral, coronary,
pulmonary). From left to right: $\mathcal{R}$ (green), $\mathcal{L}$ (orange),
$\mathcal{C}$ (purple), and $\mathcal{S}$ (pink).

Figure 13: Correlation of optimized 0D elements from 3D data with 0D ground
truth for all $N=72$ models (aorta, aortafemoral, coronary, pulmonary). From
left to right: $\mathcal{R}$ (green), $\mathcal{L}$ (orange), $\mathcal{C}$
(purple), and $\mathcal{S}$ (pink).

## References

  * [1] E. L. Schwarz, et al.  Beyond cfd: Emerging methodologies for predictive simulation in cardiovascular health and disease.  Biophysics Reviews, 4(1):011301, 2023.
  * [2] I. E. Vignon-Clementel, et al.  Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries.  Computer Methods in Applied Mechanics and Engineering, 195(29-32):3776–3796, 2006.
  * [3] M. R. Pfaller, et al.  Automated generation of 0D and 1D reduced-order models of patient-specific blood flow.  International Journal for Numerical Methods in Biomedical Engineering, 38(10), 2022.
  * [4] P. J. Nair, et al.  Non-invasive estimation of pressure drop across aortic coarctations: Validation of 0d and 3d computational models with in vivo measurements.  Annals of Biomedical Engineering, 2024.
  * [5] J. Schrauwen, et al.  Geometry-based pressure drop prediction in mildly diseased human coronary arteries.  Journal of Biomechanics, 47(8):1810–1815, 2014.
  * [6] L. Garber, et al.  The critical role of lumped parameter models in patient-specific cardiovascular simulations.  Archives of Computational Methods in Engineering, 29(5):2977–3000, 2021.
  * [7] I. E. Vignon-Clementel, et al.  A primer on computational simulation in congenital heart disease for the clinician.  Progress in Pediatric Cardiology, 30(1-2):3–13, 2010.
  * [8] C. Chnafa, et al.  Improved reduced-order modelling of cerebrovascular flow distribution by accounting for arterial bifurcation pressure drops.  Journal of Biomechanics, 51:83–88, 2017.
  * [9] M. Mirramezani et al.  A distributed lumped parameter model of blood flow.  Annals of Biomedical Engineering, 48(12):2870–2886, 2020.
  * [10] R. Pewowaruk, et al.  Accelerated estimation of pulmonary artery stenosis pressure gradients with distributed lumped parameter modeling vs. 3d CFD with instantaneous adaptive mesh refinement: Experimental validation in swine.  Annals of Biomedical Engineering, 49(9):2365–2376, 2021.
  * [11] I. S. Lan, et al.  Virtual transcatheter interventions for peripheral pulmonary artery stenosis in williams and alagille syndromes.  Journal of the American Heart Association, 11(6), 2022.
  * [12] J. D. Lee, et al.  A probabilistic neural twin for treatment planning in peripheral pulmonary artery stenosis.  arXiv, 2023.
  * [13] J. S. Tran, et al.  Automated tuning for parameter identification and uncertainty quantification in multi-scale coronary simulations.  Computers & Fluids, 142:128–138, 2017.
  * [14] A. L. Marsden.  Optimization in cardiovascular modeling.  Annual Review of Fluid Mechanics, 46(1):519–546, 2014.
  * [15] G. Ninos, et al.  Uncertainty quantification implementations in human hemodynamic flows.  Computer Methods and Programs in Biomedicine, 203:106021, 2021.
  * [16] J. Nitzler, et al.  A generalized probabilistic learning approach for multi-fidelity uncertainty quantification in complex physical simulations.  Computer Methods in Applied Mechanics and Engineering, 400:115600, 2022.
  * [17] J. Biehler, et al.  Towards efficient uncertainty quantification in complex and large-scale biomechanical problems based on a bayesian multi-fidelity scheme.  Biomechanics and modeling in mechanobiology, 14:489–513, 2015.
  * [18] J. Biehler, et al.  Multifidelity approaches for uncertainty quantification.  GAMM-Mitteilungen, 42(2):e201900008, 2019.
  * [19] P.-S. Koutsourelakis.  Accurate uncertainty quantification using inaccurate computational models.  SIAM Journal on Scientific Computing, 31(5):3274–3300, 2009.
  * [20] D. Nolte et al.  Inverse problems in blood flow modeling: A review.  International Journal for Numerical Methods in Biomedical Engineering, 38(8), 2022.
  * [21] R. L. Spilker et al.  Tuning multidomain hemodynamic simulations to match physiological measurements.  Annals of Biomedical Engineering, 38(8):2635–2648, 2010.
  * [22] M. Ismail, et al.  Adjoint-based inverse analysis of windkessel parameters for patient-specific vascular models.  Journal of Computational Physics, 244:113–130, 2013.
  * [23] K. Menon, et al.  Predictors of Myocardial Ischemia in Patients with Kawasaki Disease: Insights from Patient-Specific Simulations of Coronary Hemodynamics.  Journal of Cardiovascular Translational Research, 16(5):1099–1109, 2023.
  * [24] K. Menon, et al.  Personalized coronary and myocardial blood flow models incorporating ct perfusion imaging and synthetic vascular trees.  medRxiv, 2023.
  * [25] S. Pant, et al.  A methodological paradigm for patient-specific multi-scale CFD simulations: From clinical measurements to parameter estimates for individual analysis.  International Journal for Numerical Methods in Biomedical Engineering, 30(12):1614–1648, 2014.
  * [26] M. Spitieris, et al.  Bayesian calibration of Arterial Windkessel Model, 2022.  arXiv:2201.06883 [stat].
  * [27] J. Seo, et al.  Multifidelity estimators for coronary circulation models under clinically informed data uncertainty.  International Journal for Uncertainty Quantification, 10(5):449–466, 2020.
  * [28] C. M. Fleeter, et al.  Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics.  Computer Methods in Applied Mechanics and Engineering, 365:113030, 2020.
  * [29] D. P. Bertsekas et al.  Introduction to probability.  Athena Scientific, 2002.
  * [30] J. Kaipio et al.  Statistical and computational inverse problems.  Number v. 160 in Applied mathematical sciences. Springer, New York, 2005\.
  * [31] N. Chopin et al.  An Introduction to Sequential Monte Carlo.  Springer, 2020.
  * [32] N. Chopin.  A sequential particle filter method for static models.  Biometrika, 89(3):539–552, 2002.  Publication Title: Biometrika Volume: 89 Issue: 3.
  * [33] I. E. Vignon-Clementel, et al.  Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries.  Computer Methods in Applied Mechanics and Engineering, 195(29-32):3776–3796, 2006.
  * [34] I. Vignon-Clementel, et al.  Outflow boundary conditions for 3d simulations of non-periodic blood flow and pressure fields in deformable arteries.  Computer Methods in Biomechanics and Biomedical Engineering, 13(5):625–640, 2010.
  * [35] M. E. Moghadam, et al.  A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations.  Journal of Computational Physics, 244:63–79, 2013.
  * [36] M. R. Pfaller, et al.  On the periodicity of cardiovascular fluid dynamics simulations.  Annals of Biomedical Engineering, 49(12):3574–3592, 2021.
  * [37] J. P. Mynard et al.  A unified method for estimating pressure losses at vascular junctions.  International Journal for Numerical Methods in Biomedical Engineering, 31(7), 2015.
  * [38] N. L. Rubio, et al.  Hybrid physics-based and data-driven modeling of vascular bifurcation pressure differences.  arXiv, 2024.
  * [39] svSolver (https://github.com/SimVascular/svSolver), 2024.
  * [40] M. Esmaily-Moghadam, et al.  A new preconditioning technique for implicitly coupled multidomain simulations with applications to hemodynamics.  Computational Mechanics, 52(5):1141–1152, 2013.
  * [41] svZeroDSolver (https://github.com/SimVascular/svZeroDSolver), 2024.
  * [42] svSuperEstimator (https://github.com/StanfordCBCL/svSuperEstimator), 2024\.
  * [43] particles (https://github.com/nchopin/particles), 2023.
  * [44] Vascular Model Repository (https://vascularmodel.com), 2024.
  * [45] L. Pegolotti, et al.  Learning reduced-order models for cardiovascular simulations with graph neural networks.  Computers in Biology and Medicine, 168:107676, 2024.
  * [46] P. J. Nair, et al.  Hemodynamics in patients with aortic coarctation: A comparison of in vivo 4d-flow MRI and FSI simulation.  In Functional Imaging and Modeling of the Heart, pages 515–523. Springer Nature Switzerland, 2023.
  * [47] K. Bäumler, et al.  Fluid–structure interaction simulations of patient-specific aortic dissection.  Biomechanics and Modeling in Mechanobiology, 19(5):1607–1628, 2020\.
  * [48] A. Boccadifuoco, et al.  Impact of uncertainties in outflow boundary conditions on the predictions of hemodynamic simulations of ascending thoracic aortic aneurysms.  Computers & Fluids, 165:96–115, 2018.
  * [49] M. N. Antonuccio, et al.  Effects of uncertainty of outlet boundary conditions in a patient-specific case of aortic coarctation.  Annals of Biomedical Engineering, 49(12):3494–3507, 2021.
  * [50] J. Chung et al.  A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-$\alpha$ method.  Journal of Applied Mechanics, 60(2):371–375, 1993.
  * [51] K. E. Jansen, et al.  A generalized-alpha method for integrating the filtered navier–stokes equations with a stabilized finite element method.  Computer Methods in Applied Mechanics and Engineering, 190(3-4):305–319, 2000.
  * [52] K. Levenberg.  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics, 2(2):164–168, 1944.
  * [53] D. W. Marquardt.  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the Society for Industrial and Applied Mathematics, 11(2):431–441, 1963.