# Physiologically valid 3D-0D closed loop model of the heart and circulation -
Modeling the acute response to altered loading and contractility

Christoph M. Augustin  Matthias A. F. Gsell  Elias Karabelas  Gernot Plank
Gottfried Schatz Research Center: Division of Biophysics, Medical University
of Graz, Graz, Austria  BioTechMed-Graz, Graz, Austria

###### Abstract

Computer models of cardiac electro-mechanics (EM) have the potential to
enhance our understanding of underlying mechanisms in the analysis of normal
and dysfunctional hearts and show promise as a guidance in designing and
developing new treatment. Here, pre- and afterload conditions play a critical
role in governing and regulating heart function under both normal and
pathological conditions. The common approach of imposing loads upon the heart
by isolated pre- and afterload models does not fully address all properties of
the circulatory system as important feedback mechanisms are missing. Closed-
loop models are capable of modeling these features and also guarantee the
conservation of blood volume throughout the circulatory system over multiple
cardiac cycles.

In this study, we present the coupling of a 3D bi-ventricular EM model to the
sophisticated 0D closed-loop model CircAdapt. We show the capabilities of this
advanced framework by comparing a control case of a healthy patient to
different pathological scenarios. These simulations setups are not feasible
using isolated pre- and afterload or simple closed-loop models and are
hitherto unprecedented in the literature.

###### keywords:

3D-0D coupling, cardiovascular biomechanics, patient-specific modeling cardiac
electromechanics.

##  1 Introduction

Cardiovascular diseases (CVDs) are the primary cause of mortality and
morbidity in industrialized nations and pose a significant burden on health
care systems worldwide [1, 2, 3]. Despite continuous advances in diagnostic
and therapeutic modalities that provide a wealth of medical data, the optimal
treatment of CVDs remains a major challenge [4]. In no small part, this is due
to the highly complex multiphysics nature of cardiac function — the heart is
an electrically controlled mechanical pump driving blood through the
circulatory system — but also due to the lack of efficient tools which would
allow the comprehensive analysis of large and disparate datasets acquired in
the clinic. Computer models of cardiac electro-mechanics (EM) show high
promise to fill this void. Computational modeling and simulation of cardiac EM
offers a framework that is able to capture mechanistic relations between
clinical observations quantitatively. As such, they are increasingly used as a
clinical research tool and, potentially, will evolve in the future towards a
clinical modality that aides in diagnosis, therapy stratification and planning
to better inform clinical decision making. In particular, mechanistic patient-
specific cardiac EM model may be able to stratify therapies and make
predictions on the acute outcomes of therapies altering EM function such as
valve repair/replacement, stenting procedures, or pacing therapies such as
CRT. Their ability to predict acute therapeutic outcomes, once proven, will
pave the way towards precision cardiology where therapies are tailored for
individual patients .

Such advanced clinical applications of cardiac EM models critically depend on
their ability i) to comprehensively represent mechanisms underlying a broader
range of physiological EM function, and, ii) to tailor models to approximate –
with acceptable fidelity – anatomy and physiology in a patient-specific
manner. Such comprehensive EM models are of significant complexity as the
interplay between all major components governing a heart beat must be taken
into account. Key components comprise a model of cardiac electophysiology (EP)
to produce electrical activation and repolarization patterns, cardiac
mechanics to describe deformation, and the coupling between these physics by
myofilament models accounting for active contraction and relaxation under
given loading conditions imposed by a circulatory system.

The circulatory system as an extracardiac factor imposes a pressure and volume
load upon the heart and, thus, influences its mechanical activity. The
interaction between circulatory system and the heart is bidirectional, that
is, flow of blood in and out of a cavity and pressure in the cavity depend on
the current state of the circulatory system and pressure and flow in the
cardiovascular system are determined by the current state of all cavities..
The full physics of this interaction is most accurately posed as a fluid-
structure interaction problem where pressure and flow velocity of the fluid
are the coupling variables [5, 6]. Such distributed PDE-based FSI models are
most relevant for investigating local physical effects such as the flow
patterns due to deformation. However, at the level of the whole circulatory
physiology these are less suitable. As the more global pressure and flow
behavior within the cardiovascular system is governing overall EM function of
the heart, simpler lumped models have been preferred to provide appropriate
loading conditions to the heart.

Ventricular afterload is typically represented by lumped 0D Windkessel-type
models comprising 2-, 3- or 4-elements [7, 8, 9, 10, 11, 12, 13]. Less common
are more advanced 1D models derived from Navier–Stokes equations [14, 15, 16,
17, 18]. These account for effects due to pulse wave transmission and
reflection, but identifying and constraining model parameters based on
invasive measurement of the blood pressure and image-based flow measurements
is substantially more challenging than for 0D models [19].

However, while such isolated models of ventricular pre- and afterload may well
approximate the immediate acute response for a single beat, they are less well
suited for modeling transient behaviors in multi-beat simulations. Due to the
lack of coupling inconsistencies inevitably arise, leading to discrepancies in
stroke volumes between the chambers. Due to systemic coupling any alteration
of load or active contractile performance in one chamber causes a fluid shift
between compartments of the circulatory system. These would balanced out
transient differences in stroke volume until the overall cardiovascular system
settles in a new limit cycle. As an example, pulmonary hypertension associated
with left heart (mitral valve) disease may lead to a pressure-overloaded right
ventricle, ultimately causing right heart failure [20, 21]. Closed-loop
systems [22, 23, 24, 25, 26] are capable of modeling these important feedback
mechanisms and ascertain the conservation of blood volume throughout the
circulatory system over multiple cardiac cycles.

On physiological grounds therefore, closed-loop circulatory models are needed
to capture transient hemodynamic mechanisms. However, the robust and efficient
coupling of a 3D non-linear PDE-based model of cardiac EM to a lumped 0D
closed-loop model of the circulatory system remains to be a challenging
problem. The lumped hydrostatic pressure in the cavity and the flow between
cavity and circulatory system serve as coupling variables that act as pressure
boundary condition and impose volume constraints on the PDE-based cavity
models. Previous studies addressed 3D solid-0D fluid coupled problems using a
simple partitioned [27, 28, 29] or a more advanced strongly coupled monolithic
approach [30, 31, 32, 33, 19, 34]. Yet, there are only a few reports on the
coupling of a closed-loop 0D model to a 3D solid model and these consider
mostly the coupling to a relatively simple electrical equivalent circuit model
[32, 28, 33], but not to more advanced 0D state-of-the-art framework such as
the non-linear lumped parameter model CircAdapt [22] which models circulation
dynamics by a smart combination of known physiological principles.

Based on our previous work on anatomically accurate EM models of the heart
[35, 36, 37] we report on the development of a monolithic 3D solid - 0D fluid
coupling of a 3D bi-ventricular PDE model to an advanced widely used cicrAdapt
0D closed-loop model [38]. A detailed description of the underlying numerical
methods is given, including a detailed and complete mathematical description
of the circAdapt model that has been lacking so far. Using a comprehensive set
of experimental data we build a detailed canine biventricular EM model and
couple with circAdapt representing the atrial chambers and the remaining
circulatory system. After parameterization of the coupled model under baseline
conditions we first investigate efficiency and robustness of numerical scheme
and solver implementation. To rigorously validate the model its ability to
correctly replicate known physiological behaviors under transient alterations
of loading conditions and contractility is demonstrated. Finally, the ability
to replicate the standard multibeat detection protocols used for the
experimental determination of contractility is demonstrated.

##  2 Methodology

###  2.1 Biventricular FE Models

In a previous study mongrel dog data was acquired to show the maintenance
cardiac contractile coordination, pump function and efficiency of LV septal
and LV apical pacing, see [39]. The animals were handled according to the
Dutch Law on Animal Experimentation (WOD) and the European Directive for the
Protection of Vertebrate Animals Used for Experimental and Other Scientific
Purposes (86/609/EU). The protocol was approved by the Maastricht University
Experimental Animal Committee.

Images were acquired on a Philips Gyroscan $1.5\text{\,}\mathrm{T}$ (NT,
Philips Medical Systems, Best, the Netherlands). The RF receiver coil was a
standard synergy body coil for thorax examinations. Images of seven short-axis
cross-sections, slice thickness $8\text{\,}\mathrm{mm}$ with inter-slice
distance $0\text{\,}\mathrm{mm}$, were obtained to capture the whole heart. LV
pressure and volume were determined using the conductance catheter technique
(CD-Leycom, The Netherlands), see [40], and the signals were digitized at
$1\text{\,}\mathrm{kHz}$.

Multilabel BiV segmentations were generated from the seven MRI cross sections
using Seg3D [41]. Each cross section was first semi-automatically segmented
using thresholding techniques with manual correction. A coherent and smooth
BiV segmentation was then generated by using an automated iterative erosion
and dilation smoothing scheme.

Multilabel meshes were created from the segmentations using the Computational
Geometry Algorithms Library, CGAL (www.cgal.org) and subsequently smoothed
with Meshtool [42]. A rule-based method according to [43] was applied to
define fiber and sheet architecture, with fiber angles changing linearly from
$-60^{\circ}$ at the epicardium to $+60^{\circ}$ at the endocardium
[streeter]. Universal ventricular coordinates were computed [bayer2018:_uvc]
to support the flexible definition of stimulation sites and mechancial
boundary conditions. For a single model a coarse mesh was generated to reduce
the computational cost. Average edge lengths were chosen to ascertain that at
least two elements were generated transmurally across the myocardial walls.
The average edge length in the LV is $\sim\\!$3.4\text{\,}\mathrm{mm}$$ and in
the RV $\sim\\!$2.4\text{\,}\mathrm{mm}$$, see Fig. 1.

Figure 1: Model discretization and fibers assigned to the model.

###  2.2 Electromechanical PDE Model

#### Tissue Mechanics

Cardiac tissue is mechanically characterized as a hyperelastic, nearly
incompressible, orthotropic material with a nonlinear stress-strain
relationship. The deformation gradient $\bm{\mathbf{F}}$ describes the
deformation ${\mathbf{u}}$ of a body from the reference configuration
$\Omega_{0}({\mathbf{X}})$ to the current configuration
$\Omega_{t}({\mathbf{x}})$,

$F_{ij}=\frac{\partial x_{i}}{\partial X_{j}},\quad i,j=1,2,3.$ (1)

By convention, we denote $J=\det\bm{\mathbf{F}}>0$ and introduce the right
Cauchy–Green tensor $\bm{\mathbf{C}}=\bm{\mathbf{F}}^{\top}\bm{\mathbf{F}}$.
The nearly incompressible behavior is modeled by a multiplicative
decomposition of the deformation gradient, see, e.g., [44], of the form

$\bm{\mathbf{F}}=J^{1/3}\overline{\bm{\mathbf{F}}},\
\bm{\mathbf{C}}=J^{2/3}\overline{\bm{\mathbf{C}}},\ \text{with
}\det\overline{\bm{\mathbf{F}}}=\det\overline{\bm{\mathbf{C}}}=1.$ (2)

Mechanical deformation is described by the equilibrium equations given as

$\rho_{0}\ddot{{\mathbf{u}}}(t,{\mathbf{X}})-\operatorname{Div}\left[\bm{\mathbf{F}}\bm{\mathbf{S}}({\mathbf{u}},{\mathbf{X}})\right]={\mathbf{0}}\quad\mbox{for}\;{\mathbf{X}}\in\Omega_{0}\times(0,T),$
(3)

with initial conditions

$\displaystyle{\mathbf{u}}({\mathbf{X}},0)={\mathbf{0}},\quad\dot{{\mathbf{u}}}({\mathbf{X}},0)={\mathbf{0}}.$

Here, $\rho_{0}$ is the density in reference configuration;
$\ddot{{\mathbf{u}}}$ are nodal accelerations; $\dot{{\mathbf{u}}}$ are nodal
velocities; $\bm{\mathbf{S}}({\mathbf{u}},{\mathbf{X}})$ is the second
Piola–Kirchhoff stress tensor; and $\operatorname{Div}$ denotes the divergence
operator in the reference configuration.

The boundary of the bi-ventricular models was decomposed in three parts,
$\partial\Omega_{0}=\overline{\Gamma}_{\mathrm{endo},0}\cup\overline{\Gamma}_{\mathrm{epi},0}\cup\overline{\Gamma}_{\mathrm{base},0}$,
with $\overline{\Gamma}_{\mathrm{endo},0}$ the endocardium,
$\overline{\Gamma}_{\mathrm{epi},0}$ the epicardium, and
$\overline{\Gamma}_{\mathrm{base},0}$ the base of the ventricles.

Normal stress boundary conditions were imposed on the endocardium

$\bm{\mathbf{F}}\bm{\mathbf{S}}({\mathbf{u}},{\mathbf{X}})\,{\mathbf{n}^{\mathrm{out}}_{0}}({\mathbf{X}})=-p(t)J\bm{\mathbf{F}}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}({\mathbf{X}})\quad\text{on}\quad{\Gamma}_{\mathrm{endo},0}\times\times(0,T)$
(4)

with ${\mathbf{n}^{\mathrm{out}}_{0}}$ the outer normal vector; omni-
directional spring type boundary conditions constrained the ventricles at the
basal cut plane $\overline{\Gamma}_{\mathrm{base},0}$ [45]; and to simulate
the mechanical constrains imposed by the pericardium spatially varying normal
Robin boundary conditions were applied at the epicardium
$\overline{\Gamma}_{\mathrm{epi},0}$ [46].

Apart from external loads the deformation of cardiac tissue is in particular
governed by active stresses intrinsically generated during contraction. To
simulate both the active and passive properties of the tissue, the total
stress $\bm{\mathbf{S}}$ is additively decomposed according to

$\bm{\mathbf{S}}=\bm{\mathbf{S}}_{\mathrm{p}}+\bm{\mathbf{S}}_{\mathrm{a}},$
(5)

where $\bm{\mathbf{S}}_{\mathrm{p}}$ and $\bm{\mathbf{S}}_{\mathrm{a}}$ refer
to the passive and active stresses, respectively.

#### Passive Stress

Passive stresses are modeled based on the constitutive equation

$\bm{\mathbf{S}}_{p}=2\frac{\partial\Psi(\bm{\mathbf{C}})}{\partial\bm{\mathbf{C}}},$
(6)

where $\Psi$ is an invariant-based strain-energy function to model the
orthotropic behavior of cardiac tissue. The prevailing orientation of
myocytes, referred to as fiber orientation, is denoted as ${\mathbf{f}}_{0}$.
Individual myocytes are surrounded and interconnected by collagen, forming
sheets, which is described by the sheet orientation ${\mathbf{s}}_{0}$,
perpendicular to ${\mathbf{f}}_{0}$. Together with the sheet-normal axis
${\mathbf{n}}_{0}$, orthogonal to the sheet and the fiber orientations, this
forms a right-handed orthonormal set of basis vectors.

Following Usyk et al. [47] the orthotropic constitutive relation is defined as

$\Psi(\bm{\mathbf{C}})=\frac{\kappa}{2}{\left(\log\,J\right)}^{2}+\frac{a}{2}\left[\exp(\mathcal{Q})-1\right],$
(7)

where the first term is the volumetric energy with the bulk modulus
$\kappa\gg$0\text{\,}\mathrm{kPa}$$ which penalizes local volume changes to
enforce near incompressible behavior of the tissue and the term in the
exponent is

$\mathcal{Q}=b_{\mathrm{ff}}\overline{E}_{\mathrm{ff}}^{2}+b_{\mathrm{ss}}\overline{E}_{\mathrm{ss}}^{2}+b_{\mathrm{nn}}\overline{E}_{\mathrm{nn}}^{2}+b_{\mathrm{fs}}\left(\overline{E}_{\mathrm{fs}}^{2}+\overline{E}_{\mathrm{sf}}^{2}\right)+b_{\mathrm{fn}}\left(\overline{E}_{\mathrm{fn}}^{2}+\overline{E}_{\mathrm{nf}}^{2}\right)+b_{\mathrm{ns}}\left(\overline{E}_{\mathrm{ns}}^{2}+\overline{E}_{\mathrm{sn}}^{2}\right).$
(8)

Here, the directional strains read

$\displaystyle\overline{E}_{\mathrm{ff}}$
$\displaystyle={\mathbf{f}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{f}}_{0},\
\overline{E}_{\mathrm{ss}}={\mathbf{s}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{s}}_{0},\
\overline{E}_{\mathrm{nn}}={\mathbf{n}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{n}}_{0},\
\overline{E}_{\mathrm{fs}}={\mathbf{f}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{s}}_{0},\
\overline{E}_{\mathrm{fn}}={\mathbf{f}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{n}}_{0},$
$\displaystyle\overline{E}_{\mathrm{ns}}$
$\displaystyle={\mathbf{n}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{s}}_{0},\
\overline{E}_{\mathrm{sf}}={\mathbf{s}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{f}}_{0},\
\overline{E}_{\mathrm{nf}}={\mathbf{n}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{f}}_{0},\
\overline{E}_{\mathrm{sn}}={\mathbf{s}}_{0}\cdot\overline{\bm{\mathbf{E}}}{\mathbf{n}}_{0},$

with
$\overline{\bm{\mathbf{E}}}=\frac{1}{2}(\overline{\bm{\mathbf{C}}}-\bm{\mathbf{I}})$
the modified isochoric Green–Lagrange strain tensor and parameters given in
Tab. 1.

#### Active Stress

Stresses due to active contraction are assumed to be orthotropic with full
contractile force along the myocyte fiber orientation ${\mathbf{f}}_{0}$ and
$40\text{\,}\%$ contractile force along the sheet orientation
${\mathbf{s}}_{0}$ [48, 49]. Thus, the active stress tensor is defined as

$\bm{\mathbf{S}}_{\mathrm{a}}=S_{\mathrm{a}}{\left({\mathbf{f}}_{0}\cdot\bm{\mathbf{C}}{\mathbf{f}}_{0}\right)}^{-1}{\mathbf{f}}_{0}\otimes{\mathbf{f}}_{0}+$0.4$\,S_{\mathrm{a}}{\left({\mathbf{s}}_{0}\cdot\bm{\mathbf{C}}{\mathbf{s}}_{0}\right)}^{-1}{\mathbf{s}}_{0}\otimes{\mathbf{s}}_{0},$
(9)

where $S_{\mathrm{a}}$ is the scalar active stress describing the contractile
force. A simplified phenomenological contractile model was used to represent
active stress generation [50]. Owing to its small number of parameters and its
direct relation to clinically measurable quantities such as peak pressure, and
the maximum rate of rise of pressure this model is fairly easy to fit and thus
very suitable for being used in clinical EM modeling studies. Briefly, the
active stress transient is given by

$S_{\mathrm{a}}(t,\lambda)=S_{\mathrm{peak}}\,\phi(\lambda)\,\tanh^{2}\left(\frac{t_{\mathrm{s}}}{\tau_{\mathrm{c}}}\right)\,\tanh^{2}\left(\frac{t_{\mathrm{dur}}-t_{\mathrm{s}}}{\tau_{\mathrm{r}}}\right),\qquad\text{for
}0<t_{\mathrm{s}}<t_{\mathrm{dur}},$ (10)

with

$\phi=\tanh(\mathrm{ld}(\lambda-\lambda_{0})),\quad\tau_{\mathrm{c}}=\tau_{\mathrm{c_{0}}}+\mathrm{ld}_{\mathrm{up}}(1-\phi),\quad
t_{\mathrm{s}}=t-t_{\mathrm{a}}-t_{\mathrm{emd}}$ (11)

and $t_{\mathrm{s}}$ is the onset of contraction; $\phi(\lambda)$ is a non-
linear length-dependent function in which $\lambda$ is the fiber stretch and
$\lambda_{0}$ is the lower limit of fiber stretch below which no further
active tension is generated; $t_{\mathrm{a}}$ is the local activation time
from Eq. 12, defined when the local transmembrane potential passes the
threshold voltage $V_{\mathrm{m,thresh}}$; $t_{\mathrm{emd}}$ is the EM delay
between the onsets of electrical depolarization and active stress generation;
$S_{\mathrm{peak}}$ is the peak isometric tension; $t_{\mathrm{dur}}$ is the
duration of active stress transient; $\tau_{\mathrm{c}}$ is time constant of
contraction; $\tau_{\mathrm{c_{0}}}$ is the baseline time constant of
contraction; $\mathrm{ld}_{\mathrm{up}}$ is the length-dependence of
$\tau_{\mathrm{c}}$; $\tau_{\mathrm{r}}$ is the time constant of relaxation;
and $\mathrm{ld}$ is the degree of length dependence. For the parameter values
used in the simulations see Tab. 1. Note that active stresses in this
simplified model are only length-dependent, but dependence on fiber velocity,
$\dot{\lambda}$, is ignored.

#### Electrophysiology

A recently developed reaction-eikonal (R-E) model [51] was employed to
generate electrical activation sequences which serve as a trigger for active
stress generation in cardiac tissue. The hybrid R-E model combines a standard
reaction-diffusion (R-D) model based on the monodomain equation with an
eikonal model. Briefly, the eikonal equation is given as

$\left\\{\begin{array}[]{rcll}\sqrt{\nabla_{{\mathbf{X}}}t_{\mathrm{a}}^{\top}\,\bm{\mathbf{V}}\,\nabla_{{\mathbf{X}}}t_{\mathrm{a}}}&=&1&\text{in
}\Omega_{0},\\\ t_{\mathrm{a}}&=&t_{0}&\text{on
}\Gamma_{0}^{\ast},\end{array}\right.$ (12)

where $(\nabla_{{\mathbf{X}}})$ is the gradient with respect to the end-
diastolic reference configuration $\Omega_{0}$; $t_{\mathrm{a}}$ is a positive
function describing the wavefront arrival time at location
${\mathbf{X}}\in\Omega_{0}$; and $t_{0}$ are initial activations at locations
$\Gamma_{0}^{\ast}\subseteq\Gamma_{\mathrm{N},0}$. The symmetric positive
definite $3\times 3$ tensor $\bm{\mathbf{V}}({\mathbf{X}})$ holds the squared
velocities
$\left(v_{\mathrm{f}}({\mathbf{X}}),v_{\mathrm{s}}({\mathbf{X}}),v_{\mathrm{n}}({\mathbf{X}})\right)$
associated to the tissue’s eigenaxes ${\mathbf{f}}_{0}$, ${\mathbf{s}}_{0}$,
and ${\mathbf{n}}_{0}$. The arrival time function
$t_{\mathrm{a}}({\mathbf{X}})$ was subsequently used in a modified monodomain
R-D model given as

$\beta C_{\mathrm{m}}\frac{\partial V_{\mathrm{m}}}{\partial
t}=\nabla_{{\mathbf{X}}}\cdot\bm{\sigma}_{\mathrm{i}}\nabla_{{\mathbf{X}}}V_{\mathrm{m}}-\beta
I_{\mathrm{ion}}+I_{\mathrm{foot}},$ (13)

with $\beta$ the membrane surface-to-volume ratio; $C_{\mathrm{m}}$ the
membrane capacitance; $V_{\mathrm{m}}$ the unknown transmembrane voltage;
$\bm{\sigma}_{\mathrm{i}}$ the intracellular conductivity tensor which holds
the scalar conductivities ($g_{\mathrm{f}}(\bm{X}),\ g_{\mathrm{s}}(\bm{X}),\
g_{\mathrm{n}}(\bm{X})$) and is coupled to $\bm{\mathbf{V}}({\mathbf{X}})$
proportionally [52]; and $I_{\mathrm{ion}}$ the membrane ionic current
density. Additionally, an arrival time dependent foot current,
$I_{\mathrm{foot}}(t_{\mathrm{a}})$, was added which is designed to mimic
subthreshold electrotonic currents to produce a physiological foot of the
action potential. The key advantage of the R-E model is its ability to compute
activation sequences at much coarser spatial resolutions that are not
afflicted by the spatial undersampling artifacts leading to conduction slowing
or even numerical conduction block as it is observed in standard R-D models.
Ventricular EP was represented by the tenTusscher–Noble–Noble–Panfilov model
of the human ventricular myocyte [53].

#### Computation of Volumes

To compute the flow into the cardiovascular system, the volume of each cavity
— i.e., left ventricle (lv), left atrium (la), right ventricle (rv), right
atrium (ra) — that is modeled as a 3D elastic body has to be tracked as a
function of time: $V^{\mathrm{PDE}}_{c}({\mathbf{x}},t)$ for
$c\in\\{\mathrm{lv,rv,la,ra}\\}$. A reduction in cavitary volume

$\frac{\partial V^{\mathrm{PDE}}_{c}({\mathbf{x}},t)}{\partial t}<0,$

drives a positive flow into the cavitary system. In a purely electro-
mechanical simulation context where the fluid domain is not modeled
explicitly, the cavitary volume is not discretized, only the surface
$\Gamma_{c}$ enclosing the cavitary volume is known. Assuming that the entire
surface of the cavitary volume is available, that is, also the faces
representing the valves are explicitly discretized, the enclosed volume
$V^{\mathrm{PDE}}_{c}$ can be computed from this surface using the divergence
theorem

$V^{\mathrm{PDE}}_{c}({\mathbf{u}},t)=V^{\mathrm{PDE}}_{c}({\mathbf{x}},t)=\frac{1}{3}\int_{\Gamma_{c,t}}{\mathbf{x}}\cdot{\mathbf{n}}\,\mathrm{d}\Gamma_{c,t}.$
(14)

Using this approach, the volume $V^{\mathrm{PDE}}_{c}({\mathbf{x}},t)$ can be
computed for each state of deformation at time $t$ and the flow can be derived
by differentiating with respect to $t$.

###  2.3 Lumped ODE model of the circulatory system: the CircAdapt model

CircAdapt [22], as shown schematically in Fig. 2, is a lumped 0D model of the
human heart and circulation. It enables real-time simulation of cardiovascular
system dynamics in a wide variety of physiological and pathophysiological
situations. The entire cardiovascular system is modeled as a concatenation of
modules: a tube module representing the systemic and pulmonary arteries and
veins (A.1); a chamber module modeling actively contracting chambers, i.e.,
left and right atria and left and right ventricles (A.3) where myofiber
mechanics and contraction is described by a sarcomere module (A.2); following
Lumens et al. [54] this also includes inter-ventricular mechanical interaction
through the inter-ventricular septum (A.4); a valve module representing the
aortic, mitral, pulmonary, and tricuspid valves (A.8); a module representing
systemic and pulmonary peripheral microvasculatures (A.6); and a module
accounting for effects of the pericardium (A.5). The modules are connected by
flows over valves and venous-atrial inlets (A.7) and the whole lumped model
consists of 26 ordinary differential equations (ODEs). In this work we use an
adaptive Runge–Kutta–Fehlberg method (RKF45) to solve the system of ODEs
(A.9).

In A the mathematical underpinnings of the CircAdapt model are outlined.
Briefly, cavity pressures and cavity volumes are interconnected as follows:
volumes regulate cavity wall areas, which in turn determine strain of the
myofibers in the wall. Strain is used to calculate myofiber stress, (28, 29),
which drives wall tension in each cardiac wall (41). Using Laplace’s law,
transmural pressure is calculated from wall tension and curvature for each
wall (42). Cavity pressures are found by adding the transmural pressures to
the intra-pericardial pressure surrounding the myocardial walls (50).
Consecutively, cavity pressures are used to update flow over valves (62) and
thus intra-cavitary volumes (53).

A significant advantage of the modular setup of the model is that a simple 0D
module can be straightforwardly replaced by the more complicated FE model in
Section 2.1. In this setup CircAdapt provides realistic boundary conditions to
the FE problem, see Section 2.4.

The version of the CircAdapt model used for all simulations has been published
previously [38] and can also be downloaded from the CircAdapt website
(http://www.circadapt.org).

TUBE A.1 compute $A_{t}$, $p_{t}$, $Z_{t}$ (20, 21, 22) SARCOMERE A.2 estimate
$E_{c}^{\mathrm{fib}}$, $\sigma_{c}^{\mathrm{fib}}$,
$\kappa_{c}^{\mathrm{fib}}$ (24, 30, 31) CAVITY A.3 update
$E_{c}^{\mathrm{fib}}$ (36) and $C_{c}^{\mathrm{mid}}$,
$A_{c}^{\mathrm{mid}}$, $T_{c}^{\mathrm{mid}}$ (33, 35, 41) TRISEG A.4 update
$\dot{y}^{\mathrm{mid}}$, $\dot{V}_{\mathrm{sv}}^{\mathrm{mid}}$ (43–45)
CAVITY A.3 compute $A_{c}$, $Z_{c}$, $p_{c}$ (37, 38, 42) SARCOMEREA.2 update
$\sigma_{c}^{\mathrm{fib}}$ (30)compute
$\dot{L}_{c}^{\mathrm{cont}},\dot{C}_{c}$ (26,27)  PERICARD A.5 compute
$p_{\mathrm{peri}}$ and update $p_{c}$ (49, 50) VALVEA.8compute $\dot{q}_{v}$
(62) PERIPHERYA.6 compute $q_{py}$ (51) CONNECTA.7 update
$p_{t},p_{c}$,compute $\dot{V}_{c},\dot{V}_{t}$ (52–54)  SOLVE A.9 update
$V_{t}$, $V_{c}$, $C_{c}$, $L_{c}^{\mathrm{cont}}$, $q_{v}$,
$y^{\mathrm{mid}}$, $V_{\mathrm{sv}}^{\mathrm{mid}}$ Figure 2: CircAdapt
circulatory system, based on [38], which connects tubes ($t$), cavities ($c$),
valves ($v$), and pulmonary and systemic periphery ($py$). In each timestep
the ODE system is solved using a Runge–Kutta–Fehlberg method, see A.9, to
update the ODE variables (in green), i.e., volumes of tubes ($V_{t}$) and
cavities ($V_{c}$); sarcomere contractility ($C_{c}$) and sarcomere length
($L_{c}^{\mathrm{cont}}$) for each of the cavities and the septum; flow over
valves ($q_{v}$); and septal midwall volume ($V_{\mathrm{sv}}^{\mathrm{mid}}$)
and radius ($y^{\mathrm{mid}}$), see Fig. 6c. In the following steps the
updated variables are used to compute current pressures ($p_{c},p_{t}$), cross
sectional areas ($A_{c},A_{t}$), and impedances ($Z_{c},Z_{t}$) for tubes and
cavities; fiber strain ($E_{c}^{\mathrm{fib}}$), fiber stiffness
($\kappa_{c}^{\mathrm{fib}}$), and fiber stress ($\sigma_{c}^{\mathrm{fib}}$)
for the sarcomeres of each cavity and the septum; midwall curvature
($C_{c}^{\mathrm{mid}}$), midwall area ($A_{c}^{\mathrm{mid}}$), and midwall
tension ($T_{c}^{\mathrm{mid}}$) for each cavity and the septum; pericardial
pressure $p_{\mathrm{peri}}$; and flow over the systemic and pulmonary
periphery $q_{py}$.

###  2.4 PDE-ODE Coupling

Coupling between PDE and ODE models can be achieved in various ways,
fundamentally, the problem is to find the new state of deformation
$\mathbf{u}_{n+1}$ as a function of the pressure $p_{n+1}$ in a given cavity
at time $n+1$. The pressure $p_{n+1}$ is applied as a Neumann boundary
condition at the cavitary surface. This pressure is not known and has to be
determined in a way which depends on the current state of the cavity.
Basically, two scenarios have to be considered: (i) when all valves are
closed, the cavity is in an isovolumetric state. That is, the muscle enclosing
the cavity may deform, but the volume has to remain constant. Therefore if
active stresses vary over time during an isovolumetric phase, the pressure
$p_{n+1}$ in the cavity has to vary as well to keep the cavitary volume
constant; (ii) when at least one valve is open or regurging, the cavitary
volume is changing. In this case the pressure $p_{n+1}$ is influenced by the
state of the circulatory system. Thus $p_{n+1}$ has to be determined in a way
that matches mechanical deformation and state of the system. Pressure
$p_{n+1}$ in the cardiovascular system depends on flow and flow rate which are
governed by cardiac deformation and as such the two models are tightly
bidirectionally coupled.

The simplest approach for the PDE-ODE coupling is to prescribe $p_{n+1}$
explicitly, e.g., [27, 29]. During ejection phases this is achieved by
updating cavity volumes and flow based on the current prediction on the change
in the state of deformation under the currently predicted pressure $p_{n+1}$.
In this scenario the pressure boundary condition in each non-linear solver
step is modified within each Newton iteration $k$. The new prediction
$p_{n+1}^{k+1}$ is then prescribed explicitly as a Neumann boundary condition.
While this explicit approach is easy to implement and may be incorporated into
an existing FE solver package without difficulty, it may introduce some
inaccuracies during ejection phases and it tends to be unstable during
isovolumetric phases. Instabilities stem from the problem of estimating the
change in pressure necessary to maintain the volume. Inherently, this requires
to know the pressure-volume ($p\,$-$\,V$) relation of the cavity at this given
point in time. However, this knowledge on cavitary elastance is not available
and thus iterative estimates are necessary to gradually inflate or deflate a
cavity to its prescribed volume. As the elastance properties of the cavities
are highly non-linear, an overestimation may induce oscillations and an
underestimation may lead to very slow convergence and punitively large numbers
of Newton iterations.

0D ODE system Pulmonary Veins $q_{\mathrm{po}}$ Pulmonary Outlet
$q_{\mathrm{mv}}$ Mitral Valve Aorta $q_{\mathrm{av}}$ Aortic Valve Venae
Cavae $q_{\mathrm{so}}$ Systemic Outlet $q_{\mathrm{tv}}$ Tricuspid Valve
Pulmonary Artery $q_{\mathrm{pv}}$ Pulmonary Valve
LA$p_{\mathrm{la}}$$p_{\mathrm{ao}}$$\dot{V}_{\mathrm{la}}$LVPDE$p_{\mathrm{lv}}$$\dot{V}_{\mathrm{lv}}$RA$p_{\mathrm{ra}}$$p_{\mathrm{ap}}$$\dot{V}_{\mathrm{ra}}$3D
PDE
system$p_{\mathrm{ra}}$$q_{\mathrm{in}}$$p_{\mathrm{ap}}$$q_{\mathrm{out}}$$q_{\mathrm{pv}}$$q_{\mathrm{tv}}$$p_{\mathrm{la}}$$q_{\mathrm{in}}$$p_{\mathrm{ao}}$$q_{\mathrm{out}}$$q_{\mathrm{av}}$$q_{\mathrm{mv}}$RVPDE$p_{\mathrm{rv}}$$\dot{V}_{\mathrm{rv}}$
Pulmonary Periphery Systemic Periphery Figure 3: Schematic of the electrical
equivalent circuit showing the coupling of the lumped ODE model to the 3D PDE
model. Each of the ventricles (LV, RV) and each of the atria (LA, RA) can be
modeled either as a PDE, i.e., as a 3D elastic body, or as a lumped cavity in
the CircAdapt model. In this case the ventricles are modeled as PDEs while
atria are modeled as lumped cavities. Volume changes of the cavities
$\dot{V}_{\bullet}$ are driven by flow $q_{\bullet}$ of blood over valves and
outlets. Red colors indicate oxygenated and blue colors de-oxygenated blood.

A more elaborate approach is to treat $p_{n+1}$ as an unknown [55, 28, 56]. In
addition to the equilibrium equations (3–4) this requires one further equation
for each cavity $c$ that is modeled as a 3D elastic body, see Fig. 3. Hence,
we get $N_{\mathrm{cav}}$, the number of PDE cavities, additional equations of
the form

$V_{c}^{\mathrm{PDE}}({\mathbf{u}},t)-V_{c}^{\mathrm{ODE}}(p_{c},t)=0,$ (15)

where $V_{c}^{\mathrm{PDE}}({\mathbf{u}},t)$ is the cavity volume computed as
the integral over the current surface $\Gamma_{c,t}$, see Eq. 14, and
$V_{c}^{\mathrm{ODE}}(p_{c},t)$ is the cavity volumes as predicted by the
CircAdapt model for the intra-cavitary pressure $p_{c}$, see Section 2.3 and
A.

We write
$\underline{p}_{\mathtt{c}}={[p_{c}]}_{c\in\\{\mathrm{lv,rv,la,ra}\\}}$ for
the vector of up to $1\leq N_{\mathrm{cav}}\leq 4$ pressure unknowns. Then,
linearization of the variational problem, see B.2, a Galerkin FE
discretization, see B.3, and a time integration using a generalized-$\alpha$
scheme, see C, result in solving the block system to find
$\delta\underline{u}\in\mathbb{R}^{3N}$ and
$\delta\underline{p}_{\mathtt{c}}\in\mathbb{R}^{N_{\mathrm{cav}}}$ such that

$\bm{\mathbf{K}}^{\prime}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\begin{pmatrix}\delta\underline{u}\\\
\delta\underline{p}_{\mathtt{c}}\end{pmatrix}=-\underline{K}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k}),\quad\underline{K}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k}):=\begin{pmatrix}\underline{R}_{\alpha}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\\\
\underline{R}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\end{pmatrix},$
(16)

with the updates

$\displaystyle\underline{u}^{k+1}$
$\displaystyle=\underline{u}^{k}+\delta\underline{u},$ (17)
$\displaystyle\underline{p}_{\mathtt{c}}^{k+1}$
$\displaystyle=\underline{p}_{\mathtt{c}}^{k}+\delta\underline{p}_{\mathtt{c}}.$
(18)

Here, $\underline{u}^{k}\in\mathbb{R}^{3N}$ and
$\underline{p}_{\mathtt{c}}^{k}\in\mathbb{R}^{N_{\mathrm{cav}}}$ are the
solution vectors at the $k$-th Newton step. The block tangent stiffness matrix
$\bm{\mathbf{K}}^{\prime}$ is assembled according to Eqs. 82–86 and Eq. 102
and the right hand side vector $\underline{R}_{\alpha}$ according to Eqs. 87
and 88 and Eq. 99. The residual $\underline{R}_{\mathrm{p}}$ which measures
the accuracy of the current coupling is the discrete version of (15), i.e.,

$\underline{R}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k}):=\underline{V}^{\mathrm{PDE}}(\underline{u}^{k})-\underline{V}^{\mathrm{ODE}}(\underline{p}_{\mathtt{c}}^{k}).$
(19)

The whole procedure to perform the PDE-ODE coupling is summarized in Algorithm
1.

Algorithm 1 Coupling of the lumped ODE model to the 3D PDE model

1:Initialize time $n=0$

2:Initial displacement $\underline{u}_{0}=\underline{0}$

3:Initial cavity pressures
$\underline{p}_{\mathrm{c,0}}={[p_{c}]}_{c\in\\{\mathrm{lv,rv,la,ra}\\}}$ at
time $n=0$

4:Initialize final time point $n_{\max}$ and maximal number of Newton
iterations $k_{\max}$

5:Initialize Newton tolerance $\epsilon=10^{-6}$

6:Compute initial cavity volumes
$\underline{V}^{\mathrm{PDE}}(\underline{u}^{k})={[V_{c}^{\mathrm{PDE}}(\underline{u}^{k})]}_{c\in\\{\mathrm{lv,rv,la,ra}\\}}$

7:Run CircAdapt ODE system, see Fig. 2, until steady-state is found and get
$\underline{V}^{\mathrm{ODE}}(\underline{p}_{\mathrm{c,0}})={[V_{c}^{\mathrm{ODE}}(p_{c})]}_{c\in\\{\mathrm{lv,rv,la,ra}\\}}$

8:while $n<n_{\max}$ do

9: Initialize Newton iterator: $k=0$

10: Initial guesses for Newton: $\underline{u}^{0}=\underline{u}_{n}$,
$\underline{p}_{\mathtt{c}}^{0}=\underline{p}_{\mathtt{c},{n}}$

11: while $k<k_{\max}$ do

12: Assemble block matrix $\triangleright$ Eqs. 82–86 and Eq. 102

13: and right hand side. $\triangleright$ Eqs. 87, 88, 89, 90, 99, 100 and 101

14: Solve linearized system for $\delta\underline{u}$ and
$\delta\underline{p}_{\mathtt{c}}$ $\triangleright$ Eq. 79

15: Update displacement
$\underline{u}^{k+1}=\underline{u}^{k}+\delta\underline{u}$ and cavity
pressures
$\underline{p}_{\mathtt{c}}^{k+1}=\underline{p}_{\mathtt{c}}^{k}+\delta\underline{p}_{\mathtt{c}}$

16: Update cavity volumes $\underline{V}^{\mathrm{PDE}}(\underline{u}^{k+1})$
$\triangleright$ Eq. 14

17: Update ODE system and get
$\underline{V}^{\mathrm{ODE}}(\underline{p}_{\mathtt{c}}^{k+1})$
$\triangleright$ Fig. 2

18: Convergence test:

19: if
$\lVert\underline{R}_{\alpha}(\underline{u}^{k+1},\underline{p}_{\mathtt{c}}^{k+1})\rVert_{L^{2}}<\epsilon$
and
$\lVert\underline{R}_{\mathrm{p}}(\underline{u}^{k+1},\underline{p}_{\mathtt{c}}^{k+1})\rVert_{\infty}<\epsilon$
then

20: Solution at time $n+1$: $\underline{u}_{n+1}=\underline{u}^{k+1}$,
$\underline{p}_{\mathtt{c},{n+1}}=\underline{p}_{\mathtt{c}}^{k+1}$

21: break $\triangleright$ Newton converged

22: else

23: $k=k+1$ $\triangleright$ Next Newton step

24: end if

25: end while

26: $n=n+1$$\triangleright$ Next time step

27:end while

#### Temporal synchronization of chamber contraction

In the lumped circAdapt model contraction in individual chambers is controlled
by prescribed trigger events. Based on a given heart rate contraction of the
left atrium is triggered at intervals corresponding to a basic cycle length of
$1/heartrate$. In all other chambers contraction is triggered by prescribed
delays relative to the instants of contraction of the left atrium. In a hybrid
coupled model contraction times used in 3D electromechanics (10,11) and in the
lumped CircAdapt model (27) must be synchronized accordingly. For this sake a
finite-state machine (FSM) is used to control activation cycles in both PDE
and 0D chambers. The FSM provides trigger events to the atrial cavities
implemented in the 0D model and to the stimulation units of the ventricular
cavities implemented in the PDE model. The FSM allows for a fine-grained
control of the activation sequence in the combined model, providing trigger
events for left atrium, right atrium, atrio-ventricular node entry and exit by
a given atrio-ventricular delay, left and right His bundle as well as
activation of antero-septal, septal and posterior fascicle in the LV, as well
as septal fascicle and a fascicle linked to the moderator band in the RV. The
triggers provided by the FSM can be flexibly linked within the PDE model,
either to defined stimulation events or to a His-Purkinje system to drive
electrical activation and mechanical contraction of the ventricles.

###  2.5 Numerical framework

After discretization, at each Newton–Raphson step the block system (16) has to
be solved. In that regard, we applied a Schur complement approach, see E, to
cast the problem in a pure displacement formulation and thus can reuse solver
methods established previously [35]. In brief, we used the generalized minimal
residual method (GMRES) with an relative error reduction of
$\epsilon=10^{-8}$. Efficient preconditioning was based on the library _PETSc_
[57] and the incorporated solver suite _hypre/BoomerAMG_ [58]. For the time
integration we used a generalized-$\alpha$ scheme, see C, with spectral radius
$\rho_{\infty}=0$ and damping parameters $\beta_{\mathrm{mass}}=$0.1$$,
$\beta_{\mathrm{stiff}}=$0.1$$.

We implemented the coupling scheme in the FE framework Cardiac Arrhythmia
Research Package (CARP) [59], see also the openCARP webpage
(http://www.opencarp.org). Based on the MATLAB code available on the CircAdapt
website (http://www.circadapt.org) we added a C++ circulatory system module to
CARP to guarantee a computationally efficient and strongly scaling numerical
scheme that allows fast simulation cycles.

##  3 Results

###  3.1 Model parameterization

To find a limiting cycle for the coupled 3D-0D model we started the simulation
with initializing beats of the CircAdapt model alone until a steady-state was
found for the 0D model. After that, the 3D-0D coupling scheme was activated.
To speed up simulation times we limited the number of Newton steps to
$k_{\max}=1$ for the first initial heart beats until a steady-state was found
for this semi-implicit 3D-0D model. Finally, we performed two further beats
using a fully converging Newton method with $k_{\max}=20$ and an relative
$\ell_{2}$ norm error reduction of the residual of $\epsilon=10^{-6}$.

Simulations were performed with a cardiac cycle length of
$0.585\text{\,}\mathrm{s}$ which corresponds to a heart rate (HR) of
$103\text{\,}\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$.
The pacing protocol for the FSM was as follows: RA activation time was set to
the time of the onset of the current beat; inter-atrial conduction (AA) delay
was set to $0.02\text{\,}\mathrm{s}$; atrioventricular (AV) delay to
$0.1\text{\,}\mathrm{s}$ and inter-ventricular (VV) delay to
$0\text{\,}\mathrm{s}$. Hemodynamics in the CircAdapt model reached a steady
state before it was coupled to the 3D PDE model. All parameter values used in
the CircAdapt model are listed in Tabs. 2 and 3.

For passive mechanics in the 3D PDE model initial material constants in (8)
were taken from [60] and the Klotz relation [61] was used as a target to fit
the scaling parameter $a$. Default parameters from Niederer et al. [50] for
the active contraction model (10) were manually scaled to achieve the desired
ejection fraction, peak pressure, and contraction duration. All parameters
used for the baseline case are summarized in Tab. 1.

###  3.2 Physiological validation

After parameterization of the combined EM model under baseline condition the
model’s transient response to alterations in loading conditions was
investigated to corroborate that the coupling with a closed-loop circulatory
system replicates known cardiovascular physiology. After stabilizing the model
to a baseline limit cycle step changes in preload, afterload and contractility
were imposed. In particular, left ventricular afterload was
increased/decreased by decreasing/increasing the aortic valve area by
$25\text{\,}\mathrm{\char 37\relax}$ to introduce an additional transvalvular
pressure drop. This maneuver alters arterio-ventricular coupling and the
steepness of the arterial elastance curve (4, left panel). Left ventricular
preload was then increased/decreased by decreasing/increasing the area of the
pulmonary outlet by $25\text{\,}\mathrm{\char 37\relax}$ which provokes a
transient change in venous return and thus end-diastolic filling of the
ventricle (4, right panel). Finally, ventricular contractility was
increased/decreased by increasing/decreasing the peak active stress by
$20\text{\,}\mathrm{\char 37\relax}$. This maneuver left preload and afterload
constant, but steepened the end-systolic pressure-volume relation (4, middle
panel). In addition, an instationary _walk-down_ experiment was simulated
where venous return is reduced to alter preload, as done experimentally in a
multi-beat detection protocol for estimating the load-independent cardiac
contractility indicator $E_{\rm{max}}$. The experiment was carried out for
three different states of contractility.

Figure 4: Bi-ventricular EM model coupled to CircAdapt model. Afterload is increased (decreased) by decreasing (increasing) the aortic valve area, the contractility is varied by changing the peak active stress and the preload is modified by changing the area of the pulmonary output valve. The grey area indicates ejection. Figure 5: Multibeat detection of $E_{\rm{max}}$ Parameter | Value | Unit |  Description
---|---|---|---
_Passive biomechanics_
$\rho_{0}$ | $1060.0$ | $\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{{}^{3}}$ |  tissue density
$\kappa$ | 650 | $\mathrm{kPa}$ |  bulk modulus
$a$ | 0.7 | $\mathrm{kPa}$ |  stiffness scaling
$b_{\mathrm{ff}}$ | 5.0 | [-] |  fiber strain scaling
$b_{\mathrm{ss}}$ | 6.0 | [-] |  cross-fiber in-plain strain scaling
$b_{\mathrm{nn}}$ | 3.0 | [-] |  radial strain scaling
$b_{\mathrm{fs}}$ | 10.0 | [-] |  shear strain in fiber-sheet plane scaling
$b_{\mathrm{fn}}$ | 2.0 | [-] |  shear strain in fiber-radial plane scaling
$b_{\mathrm{ns}}$ | 2.0 | [-] |  shear strain in transverse plane scaling
_Active biomechanics_
$\lambda_{0}$ | 0.7 | $\mathrm{ms}$ |  onset of contraction
$V_{\mathrm{m,Thresh}}$ | -60.0 | $\mathrm{mV}$ |  membrane potential threshold
$t_{\mathrm{emd}}$ | 15.0 | $\mathrm{ms}$ |  EM delay
$S_{\mathrm{peak}}$ | 100 (lv), 80 (rv) | $\mathrm{kPa}$ |  peak isometric tension
$t_{\mathrm{dur}}$ | 300.0 | $\mathrm{ms}$ |  duration of active contraction
$\tau_{c_{0}}$ | 100.0 | $\mathrm{ms}$ |  baseline time constant of contraction
$\mathrm{ld}$ | 5.0 | [-] |  degree of length-dependence
$\mathrm{ld}_{\mathrm{up}}$ | 500.0 | $\mathrm{ms}$ |  length-dependence of upstroke time
$\tau_{\mathrm{r}}$ | 100.0 | $\mathrm{ms}$ |  time constant of relaxation
_Electrophysiology_
$t_{\mathrm{cycle}}$ | $0.585$ | $\mathrm{s}$ |  cycle time ($=1/\mathrm{heartrate}$)
AA delay | 20.0 | $\mathrm{ms}$ |  inter-atrial conduction delay
AV delay | 100.0 | $\mathrm{ms}$ |  atrioventricular conduction delay
VV delay | 0.0 | $\mathrm{ms}$ |  inter-ventricular conduction delay
$(v_{\mathrm{f}},v_{\mathrm{s}},v_{\mathrm{n}})$ | (0.6, 0.4, 0.2) | $\mathrm{m}\mathrm{/}\mathrm{s}$ |  conduction velocities
$(v_{\mathrm{f}},v_{\mathrm{s}},v_{\mathrm{n}})$ | (1.2, 0.8, 0.4) | $\mathrm{m}\mathrm{/}\mathrm{s}$ |  fast conducting layer velocities
$(g_{\mathrm{f}},g_{\mathrm{s}},g_{\mathrm{n}})$ | (0.44, 0.54, 0.54) | $\mathrm{m}\mathrm{/}\mathrm{s}$ |  conductivities in LV and RV
$(g_{\mathrm{f}},g_{\mathrm{s}},g_{\mathrm{n}})$ | (0.65, 6.27, 6.27) | $\mathrm{m}\mathrm{/}\mathrm{s}$ |  conductivities in the fast conducting layers
$\beta$ | 1/1400 | $\mathrm{c}\mathrm{m}^{-1}$ |  membrane surface-to-volume ratio
$C_{\mathrm{m}}$ | 1 | $\mathrm{\SIUnitSymbolMicro}\mathrm{F}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$ |  membrane capacitance
Table 1: Input parameters for the 3D PDE model of the left (lv) and right (rv)
ventricle. Adjusted to match patient-specific data.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska–Curie Action
H2020-MSCA-IF-2016 InsiliCardio, GA No. 750835 and under the ERA-NET co-fund
action No. 680969 (ERA-CVD SICVALVES) funded by the Austrian Science Fund
(FWF), Grant I 4652-B to CMA. Additionally, the research was supported by the
Grants F3210-N18 and I2760-B30 from the Austrian Science Fund (FWF) and a
BioTechMed Graz flagship award “ILearnHeart” to GP.

##  Appendix A CircAdapt equations summary

###  A.1 CircAdapt Tube Module

The tube module represents the entrance of a compliant blood vessel capable of
propagating a pressure-flow wave component added to a constant flow, see [22].
Vessels directly attached to the heart, aorta (ao), arteria pulmonalis (ap),
venae cavae (vc), and venae pulmonales (vp) are modeled in a similar fashion
in CircAdapt and for the whole section we define the iterator
$t\in\\{\mathrm{ao,ap,vc,vp}\\}$.

The current lumen cross sectional area is computed by

$A_{t}=\frac{V_{t}}{l_{t}},$ (20)

with $V_{t}$ the cavity volume and $l_{t}$ the length of the vessel segment.

Using a model of a tube with a fibrous wall, see [62], gives the average
extension, $\lambda_{t}$, of the fibers in the wall by

$\lambda_{t}={\left(1+2V_{t}/V_{t}^{\mathrm{wall}}\right)}^{1/3}={\left(1+2A_{t}/A_{t}^{\mathrm{wall}}\right)}^{1/3},$

with $A_{t}^{\mathrm{wall}}$ the wall cross-sectional area and
$V_{t}^{\mathrm{wall}}=A_{t}^{\mathrm{wall}}l_{t}$ the wall volume. Cavity
pressure depends on $\lambda_{t}$

$p_{t}=\sigma_{t}(\lambda_{t})\,\lambda_{t}^{-3},$

with the mean Cauchy fiber stress $\sigma_{t}$ that is modeled by the
constitutive equation

$\sigma_{t}(\lambda_{t})=\sigma_{t}^{\mathrm{ref}}\cdot{\left(\lambda_{t}/\lambda_{t}^{\mathrm{ref}}\right)}^{k_{t}},$

see Arts et al. [62]. Here, $k_{t}$ a stiffness exponent;
$\lambda_{t}^{\mathrm{ref}}={\left(1+2V_{t}^{\mathrm{ref}}/V_{t}^{\mathrm{wall}}\right)}^{1/3}$
and
$\sigma_{t}^{\mathrm{ref}}=p_{t}^{\mathrm{ref}}\,{\left(\lambda_{t}^{\mathrm{ref}}\right)}^{3}$
are the fiber extension and fiber state at normal physiological reference
state, respectively; and $p_{t}^{\mathrm{ref}}$ is the reference tube
pressure. By combining above results the current tube pressure is computed as

$\displaystyle p_{t}$
$\displaystyle=p_{t}^{\mathrm{ref}}\,{\left(\frac{\lambda}{\lambda_{t}^{\mathrm{ref}}}\right)}^{k}\lambda_{t}^{-3}{\left(\lambda_{t}^{\mathrm{ref}}\right)}^{3}=p_{t}^{\mathrm{ref}}\,{\left(\frac{\lambda_{t}}{\lambda_{t}^{\mathrm{ref}}}\right)}^{k-3}$
$\displaystyle=p_{t}^{\mathrm{ref}}\,{\left(\frac{V_{t}^{\mathrm{wall}}+2V_{t}}{V_{t}^{\mathrm{wall}}+2V_{t}^{\mathrm{ref}}}\right)}^{\frac{k-3}{3}}=p_{t}^{\mathrm{ref}}\,{\left(\frac{A_{t}^{\mathrm{wall}}+2A_{t}}{A_{t}^{\mathrm{wall}}+2A_{t}^{\mathrm{ref}}}\right)}^{\frac{k-3}{3}},$
(21)

with $A_{t}^{\mathrm{ref}}$ the initial cross sectional area and
$V_{t}^{\mathrm{ref}}=A_{t}^{\mathrm{ref}}l_{t}$ the initial vessel volume.
The compliance is

$\displaystyle\frac{1}{C_{t}}$
$\displaystyle=\frac{\,\mathrm{d}p_{t}}{\,\mathrm{d}V_{t}}=\frac{\,\mathrm{d}}{\,\mathrm{d}V_{t}}\left[p_{t}^{\mathrm{ref}}{\left(\frac{V_{t}^{\mathrm{wall}}+2V_{t}}{V_{t}^{\mathrm{wall}}+2V_{t}^{\mathrm{ref}}}\right)}^{\frac{k-3}{3}}\right]$
$\displaystyle=p_{t}^{\mathrm{ref}}\,\frac{k-3}{3}{\left(\frac{V_{t}^{\mathrm{wall}}+2V_{t}}{V_{t}^{\mathrm{wall}}+2V_{t}^{\mathrm{ref}}}\right)}^{\left(\frac{k-3}{3}-1\right)}\frac{2}{V_{t}^{\mathrm{wall}}+2V_{t}^{\mathrm{ref}}}$
$\displaystyle=\frac{2p_{t}\,(k-3)}{3\left(V_{t}^{\mathrm{wall}}+2V_{t}\right)}=\frac{2p_{t}\,(k-3)}{3l_{t}\left(A_{t}^{\mathrm{wall}}+2A_{t}\right)}.$

Finally, the basic relation between characteristic wave impedance $Z_{t}$
compliance $C_{t}$, and inertance $I_{t}$

$Z_{t}^{2}=\frac{I_{t}}{C_{t}}=\frac{\rho_{\mathrm{b}}l_{t}}{A_{t}C_{t}}$ (22)

yields

$Z_{t}=\sqrt{\frac{2\rho_{\mathrm{b}}\,p_{t}\,l_{t}^{2}(k-3)}{3V_{t}\left(V_{t}^{\mathrm{wall}}+2V_{t}\right)}}=\sqrt{\frac{2\rho_{\mathrm{b}}\,p_{t}\,(k-3)}{3A_{t}\left(A_{t}^{\mathrm{wall}}+2A_{t}\right)}}$
(23)

with $\rho_{\mathrm{b}}$ the blood density.

###  A.2 Sarcomere mechanics

In the following a sarcomere contraction model is described that is based on a
modified Hill model, see [54, 38], for all tissue patches in the wall of the
cavity $c\in\\{\mathrm{lv,rv,sv,la,ra}\\}$. Natural strain
$E_{c}^{\mathrm{fib}}$ of the myofiber is estimated as

$E_{c}^{\mathrm{fib}}=\ln\left(\frac{L_{c}^{\mathrm{s}}}{L^{\mathrm{s,ref}}}\right)$
(24)

and from this the total sarcomere length $L_{c}^{\mathrm{s}}$ can be computed
as

$L_{c}^{\mathrm{s}}=L^{\mathrm{s,ref}}\exp\left(E_{c}^{\mathrm{fib}}\right),$
(25)

with $L_{c}^{\mathrm{s,ref}}$ a constant describing the reference sarcomere
length. The sarcomere is supposed to be made up of a contractile element of
length $L_{c}^{\mathrm{cont}}$ in series with an elastic element of length
$L_{c}^{\mathrm{elast}}=L_{c}^{\mathrm{s}}-L_{c}^{\mathrm{cont}}$.

### Sarcomere active stress

Sarcomere contracting length $L_{c}^{\mathrm{cont}}$ varies over time
according to

$\dot{L}_{c}^{\mathrm{cont}}:=\frac{\,\mathrm{d}L_{c}^{\mathrm{cont}}}{\,\mathrm{d}t}=v^{\max}\left[\frac{L_{c}^{\mathrm{s}}-L_{c}^{\mathrm{cont}}}{L^{\mathrm{elast,iso}}}-1\right],$
(26)

where the constant $L^{\mathrm{elast,iso}}$ is the length of the series
elastic element during isovolumetric contraction and the constant $v^{\max}$
is the maximum velocity of contraction.

The governing equation for the contractility $C_{c}$, describing the density
of cross bridge formation within the fibres of the patch, is

$\dot{C}_{c}:=\frac{\,\mathrm{d}C_{c}}{\,\mathrm{d}t}=f_{c}^{\mathrm{rise}}(t)C_{c}^{\mathrm{s}}\left(L_{c}^{\mathrm{cont}}\right)-f_{c}^{\mathrm{decay}}(t)C_{c}.$
(27)

In Eq. 27 sarcomere contractility rises according to the function
$f_{c}^{\mathrm{rise}}$, a phenomenological representation of the rate of
cross bridge formation within the patch,

$\displaystyle f_{c}^{\mathrm{rise}}(t)$
$\displaystyle=\frac{1}{t_{c}^{\mathrm{rise}}}0.02x_{c}^{3}{(8-x_{c})}^{2}\exp(-x_{c}),$
$\displaystyle x_{c}(t)$
$\displaystyle=\min\left(8,\max\left(0,\frac{t-t_{c}^{\mathrm{act}}}{t_{c}^{\mathrm{rise}}}\right)\right),$

depending on the time of onset of activation of the patch,
$t_{c}^{\mathrm{act}}$, and the rising time

$\displaystyle
t_{c}^{\mathrm{rise}}=0.55\tau^{\mathrm{R}}t_{c}^{\mathrm{act,ref}}.$

Here, $\tau^{\mathrm{R}}$ a constant and $t_{c}^{\mathrm{act,ref}}$ is the
reference duration of contraction for initial fiber length.

Sarcomere contractility in (27) decays according to the function
$f_{c}^{\mathrm{decay}}$

$\displaystyle f_{c}^{\mathrm{decay}}(t)$
$\displaystyle=\frac{1}{2t_{c}^{\mathrm{decay}}}\left[1+\sin\left(\operatorname{sign}(y_{c})\min\left(\frac{\pi}{2},|y_{c}|\right)\right)\right],$
$\displaystyle y_{c}(t)$
$\displaystyle=\frac{t-t_{c}^{\mathrm{act}}-t_{c}^{\mathrm{act,dur}}}{t_{c}^{\mathrm{decay}}},$

depending on the decay time

$\displaystyle
t_{c}^{\mathrm{decay}}=0.33\tau^{\mathrm{D}}t_{c}^{\mathrm{act,ref}},$

with $\tau^{\mathrm{D}}$ a constant and $t_{c}^{\mathrm{act,dur}}$ is the
duration of contraction of the fiber that lengthens with sarcomere length

$\displaystyle
t_{c}^{\mathrm{act,dur}}=\left(0.65+1.0570L_{c}^{\mathrm{norm}}\right)t_{c}^{\mathrm{act,ref}}.$

Here, $L_{c}^{\mathrm{norm}}$ is the normalized sarcomere length for active
contraction

$\displaystyle
L_{c}^{\mathrm{norm}}=\max\left(0.0001,L_{c}^{\mathrm{cont}}/L^{\mathrm{act0,ref}}-1\right),$

where $L^{\mathrm{act0,ref}}$ is the zero active stress sarcomere length.

$C_{c}^{\mathrm{s}}$ in (27) describes the increase in cross bridge formation
with intrinsic sarcomere length due to an increase in available binding sites,

$\displaystyle C_{c}^{\mathrm{s}}\left(L_{c}^{\mathrm{cont}}\right)$
$\displaystyle=\tanh\left(0.75*9.1204{\left(L_{c}^{\mathrm{norm}}\right)}^{2}\right).$

Contractility $C_{c}$ (27) and sarcomere contracting length
$L_{c}^{\mathrm{cont}}$ (26) are used to compute the actively generated fiber
stress

$\sigma_{c}^{\mathrm{fib,act}}=L^{\mathrm{act0,ref}}\sigma^{\mathrm{act,\max}}C_{c}L_{c}^{\mathrm{norm}}\frac{L_{c}^{\mathrm{s}}-L_{c}^{\mathrm{cont}}}{L^{\mathrm{elast,iso}}},$
(28)

with constants $L^{\mathrm{act0,ref}}$, $\sigma^{\mathrm{act,\max}}$,
$L^{\mathrm{elast,iso}}$, see Tab. 3.

### Sarcomere passive stress

Passive stress $\sigma_{c}^{\mathrm{fib,pas}}$ is considered to contain two
components,

$\sigma_{c}^{\mathrm{fib,pas}}=\sigma_{c}^{\mathrm{fib,tit}}+\sigma_{c}^{\mathrm{fib,ecm}},$
(29)

first the stress arising from cellular structures such as titin, a highly
abundant structural protein of the sarcomere, anchoring to the Z-disc,
$\sigma_{c}^{\mathrm{fib,tit}}$, and second the stress arising from the
extracellular matrix (ECM), $\sigma_{c}^{\mathrm{fib,ecm}}$. Both depend on
the passive fiber stretch which is computed as

$\displaystyle\lambda_{c}^{\mathrm{pas}}=\frac{L_{c}^{\mathrm{s}}}{L^{\mathrm{pas0,ref}}},$

where $L^{\mathrm{pas0,ref}}$ is sarcomere length with zero passive stress and
$L_{c}^{\mathrm{s}}$ the total sarcomere length, see above. Using that we
compute

$\displaystyle\sigma_{c}^{\mathrm{fib,tit}}=0.01\sigma^{\mathrm{act,\max}}\left({\left[\lambda_{c}^{\mathrm{pas}}\right]}^{k^{\mathrm{tit}}}-1\right),$

with $\sigma^{\mathrm{act,\max}}$ the maximal isometric stress and the
constant exponent

$\displaystyle
k^{\mathrm{tit}}=2\frac{L^{\mathrm{s,ref}}}{\,\mathrm{d}L^{\mathrm{s,pas}}}.$

The ECM is modeled as being stiffer than the myocyte contribution using

$\displaystyle\sigma_{c}^{\mathrm{fib,ecm}}=0.0349\sigma^{\mathrm{pas,\max}}\left({\left(\lambda_{c}^{\mathrm{pas}}\right)}^{10}-1\right),$

where $\sigma^{\mathrm{pas,\max}}$ is an empirical parameter.

### Sarcomere total stress

Total myofiber stress $\sigma^{\mathrm{fib}}_{c}$ is the sum of an active (28)
and a passive (29) stress component

$\sigma^{\mathrm{fib}}_{c}=\sigma_{c}^{\mathrm{fib,act}}+\sigma_{c}^{\mathrm{fib,pas}}.$
(30)

Sarcomere stiffness $\kappa_{c}^{\mathrm{fib}}$ is now computed as the
derivative of total fiber stress (30) with respect to fiber strain (24)

$\kappa_{c}^{\mathrm{fib}}=\frac{\partial\sigma_{c}^{\mathrm{fib}}}{\partial
E_{c}^{\mathrm{fib}}}=\frac{\partial\sigma_{c}^{\mathrm{fib,act}}}{\partial
E_{c}^{\mathrm{fib}}}+\frac{\partial\sigma_{c}^{\mathrm{fib,pas}}}{\partial
E_{c}^{\mathrm{fib}}},$ (31)

with

$\displaystyle\frac{\partial\sigma_{c}^{\mathrm{fib,act}}}{\partial
E_{c}^{\mathrm{fib}}}$
$\displaystyle=L^{\mathrm{act0,ref}}\sigma^{\mathrm{act,\max}}C_{c}L_{c}^{\mathrm{norm}}\frac{L_{c}^{\mathrm{s}}}{L^{\mathrm{elast,iso}}},$
$\displaystyle\frac{\partial\sigma_{c}^{\mathrm{fib,pas}}}{\partial
E_{c}^{\mathrm{fib}}}$
$\displaystyle=0.01k^{\mathrm{tit}}\sigma_{c}^{\mathrm{act,\max}}{\left(\lambda_{c}^{\mathrm{pas}}\right)}^{k^{\mathrm{tit}}}+0.0349*10\sigma^{\mathrm{pas,\max}}{\left(\lambda_{c}^{\mathrm{pas}}\right)}^{10}.$

###  A.3 CircAdapt Chamber Module

An actively contracting chamber $c\in\\{\mathrm{lv,rv,la,ra}\\}$ is modeled
using the state variables volume $V_{c}$, length of the contractile element of
the sarcomere $L_{c}^{\mathrm{cont}}$ (26), and contractility $C_{c}$ (27).
Volume changes driven by inflow and outflow of blood induce changes in midwall
volume $V_{c}^{\mathrm{mid}}$ and area $A_{c}^{\mathrm{mid}}$.

### Sphere mechanics

Note that in CircAdapt ventricles are usually modeled using the TriSeg
formulation, see A.4. If TriSeg is turned on, the calculations in this chapter
are only used for the atria while ventricular values are computed as in A.4.

Midwall volume $V_{c}^{\mathrm{mid}}$ is estimated as

$V_{c}^{\mathrm{mid}}=V_{c}+\frac{1}{2}V_{c}^{\mathrm{wall}},$ (32)

where $V_{c}^{\mathrm{wall}}$ is constant wall volume. If not set to a
specific value the wall volume is estimated by extruding the sphere enclosing
the cavity volume $V_{c}$ by a constant wall thickness
$\operatorname{h}_{c}^{\mathrm{wall}}$, see Tab. 2. Chambers are modeled as
closed spheres, thus, the following equations result from volume and surface
formulas for spheres

$\displaystyle C_{c}^{\mathrm{mid}}$
$\displaystyle={\left(\frac{4\pi}{3V_{c}^{\mathrm{mid}}}\right)}^{1/3},$ (33)
$\displaystyle A_{c}^{\mathrm{mid,tot}}$
$\displaystyle=\frac{4\pi}{{\left(C_{c}^{\mathrm{mid}}\right)}^{2}},$ (34)
$\displaystyle A_{c}^{\mathrm{mid}}$
$\displaystyle=A_{c}^{\mathrm{mid,tot}}-A_{c}^{\mathrm{mid,dead}},$ (35)

where $C_{c}^{\mathrm{mid}}$ is midwall curvature, i.e., the inverse of
radius; and $A_{c}^{\mathrm{mid,dead}}$ is non-contractile area, i.e., valve
openings and orifices.

### Update fiber strain

Natural fiber strain $E_{c}^{\mathrm{fib}}$ is calculated by

$E_{c}^{\mathrm{fib}}=\frac{1}{2}\ln\left(\frac{A_{c}^{\mathrm{mid}}}{A_{c}^{\mathrm{mid,ref}}}\right)$
(36)

with $A_{c}^{\mathrm{mid,ref}}$ the surface area in the reference state, see
[38]. Note that this updated fiber strain is used in the place of (24) to
update values in the sarcomere module A.2.

Cross-sectional area $A_{c}$ of chambers is estimated as

$\displaystyle A_{c}$
$\displaystyle=\frac{V_{c}+0.1V_{c}^{\mathrm{wall}}}{l_{c}},$ (37)
$\displaystyle l_{c}$
$\displaystyle=2{\left(V_{c}^{\mathrm{mid}}\right)}^{1/3},$

with $l_{c}$ the long-axis length of the cavity.

The characteristic wave impedance $Z_{c}$ is approximated according to (22),
see also [22], and by applying the chain rule

$Z_{c}=\frac{1}{5A_{c}}\sqrt{\rho_{\mathrm{b}}l_{c}\left|\kappa_{c}^{\mathrm{mid}}\right|},$
(38)

with the sheet stiffness

$\kappa_{c}^{\mathrm{mid}}=\frac{\partial T_{c}^{\mathrm{mid}}}{\partial
A_{c}^{\mathrm{mid}}}=\frac{V_{c}^{\mathrm{wall}}}{4{\left(A_{c}^{\mathrm{mid}}\right)}^{2}}\left(\frac{\partial\sigma_{c}^{\mathrm{fib}}}{\partial
E_{c}^{\mathrm{fib}}}-2\sigma_{c}^{\mathrm{fib}}\right)=\frac{V_{c}^{\mathrm{wall}}}{4{{\left(A_{c}^{\mathrm{mid}}\right)}^{2}}}\left(\kappa_{c}^{\mathrm{fib}}-2\sigma_{c}^{\mathrm{fib}}\right)$
(39)

and the updated fiber stiffness $\kappa_{c}^{\mathrm{fib}}$, see Eq. 31.

### Conservation of energy

CircAdapt connects midwall tension $T_{c}^{\mathrm{mid}}$ and midwall area
$A_{c}^{\mathrm{mid}}$ to fiber stress $\sigma_{c}^{\mathrm{fib}}$ and strain
$E_{c}^{\mathrm{fib}}$ through the law of conservation of energy. With the law
of Laplace we get

$\displaystyle
T_{c}^{\mathrm{mid}}\,\mathrm{d}A_{c}^{\mathrm{mid}}=\sigma_{c}^{\mathrm{fib}}V_{c}^{\mathrm{wall}}\,\mathrm{d}E_{c}^{\mathrm{fib}}$
(40)

and with (36) we get for the midwall tension

$T_{c}^{\mathrm{mid}}=\frac{\sigma_{c}^{\mathrm{fib}}V_{c}^{\mathrm{wall}}}{2A_{c}^{\mathrm{mid}}}.$
(41)

Transmural pressure $p_{c}^{\mathrm{trans}}$ is finally computed as follows

$p_{c}^{\mathrm{trans}}=2T_{c}^{\mathrm{mid}}C_{c}^{\mathrm{mid}}.$ (42)

Since at the moment external pressures are assumed to be zero, the transmural
pressure coincides with the internal pressure of the contracting chamber

$\displaystyle p_{c}=p_{c}^{\mathrm{trans}}.$

###  A.4 TriSeg model of ventricular interaction

In case that one ODE and one PDE cavity is included in the model, ventricles
are modeled as atria above. Otherwise, ventricular and septal midwall volumes
are modeled as a ventricular composite [54] which is defined by the common
radius $y^{\mathrm{mid}}$ of the wall junction and the enclosed midwall cap
volumes, see Fig. 6c. Midwall cap volumes of the right and the left ventricle
are computed as

$\displaystyle V_{\mathrm{lv}}^{\mathrm{mid}}$
$\displaystyle=-V_{\mathrm{lv}}+V^{\mathrm{mid}}_{\mathrm{sv}}-\frac{1}{2}\left(V_{\mathrm{lv}}^{\mathrm{wall}}+V_{\mathrm{sv}}^{\mathrm{wall}}\right),$
$\displaystyle V_{\mathrm{rv}}^{\mathrm{mid}}$
$\displaystyle=V_{\mathrm{rv}}+V^{\mathrm{mid}}_{\mathrm{sv}}+\frac{1}{2}\left(V_{\mathrm{rv}}^{\mathrm{wall}}+V_{\mathrm{sv}}^{\mathrm{wall}}\right).$

Here, the wall volumes of the left, $V_{\mathrm{lv}}^{\mathrm{wall}}$, and
right, $V_{\mathrm{rv}}^{\mathrm{wall}}$, ventricle are constants. The blood
pool volumes of the left, $V_{\mathrm{lv}}$, and right, $V_{\mathrm{rv}}$,
ventricle are ODE variables as well as the radius $y^{\mathrm{mid}}$ and the
septal midwall volume $V^{\mathrm{mid}}_{\mathrm{sv}}$. Note that the sign of
midwall volume $V_{c}^{\mathrm{mid}}$ is positive if wall curvature is convex
to the positive x-direction and negative otherwise.

Distance $x^{\mathrm{mid}}_{c}$, see Fig. 6c, is then computed by the relation

$V_{c}^{\mathrm{mid}}=\frac{\pi}{6}x^{\mathrm{mid}}_{c}\left({(x^{\mathrm{mid}}_{c})}^{2}+3{(y^{\mathrm{mid}})}^{2}\right),\quad\mbox{for
}c\in\\{\mathrm{lv,rv,sv}\\},$

hence

$x^{\mathrm{mid}}_{c}=q_{c}-\frac{{(y^{\mathrm{mid}})}^{2}}{q_{c}},\text{ with
}q_{c}=\sqrt[3]{\sqrt{{\left(\frac{3}{\pi}V_{c}^{\mathrm{mid}}\right)}^{2}+{\left(y^{\mathrm{mid}}\right)}^{6}}+\frac{3}{\pi}V_{c}^{\mathrm{mid}}}.$

Midwall area and curvature are consequently computed

$\displaystyle A_{c}^{\mathrm{mid}}$
$\displaystyle=\pi\left({(x^{\mathrm{mid}}_{c})}^{2}+{({y^{\mathrm{mid}}})}^{2}\right),$
$\displaystyle\mbox{for }c\in\\{\mathrm{lv,rv,sv}\\},$ $\displaystyle
C_{c}^{\mathrm{mid}}$
$\displaystyle=\frac{2x^{\mathrm{mid}}_{c}}{{(x^{\mathrm{mid}}_{c})}^{2}+{(y^{\mathrm{mid}})}^{2}},$
$\displaystyle\mbox{for }c\in\\{\mathrm{lv,rv,sv}\\},$

and used to calculate midwall tension $T_{c}^{\mathrm{mid}}$ (41).

Axial $T_{c}^{\mathrm{x}}$ and radial $T_{c}^{\mathrm{y}}$ tension components
are computed using laws of trigonometry

$\displaystyle T_{c}^{\mathrm{x}}$
$\displaystyle=T_{c}^{\mathrm{mid}}\sin\alpha,\quad\mbox{with
}\sin\alpha=\frac{2x^{\mathrm{mid}}_{c}\,y^{\mathrm{mid}}}{{(x^{\mathrm{mid}}_{c})}^{2}+{(y^{\mathrm{mid}})}^{2}},$
$\displaystyle\mbox{for }c\in\\{\mathrm{lv,rv,sv}\\},$ $\displaystyle
T_{c}^{\mathrm{y}}$
$\displaystyle=T_{c}^{\mathrm{mid}}\cos\alpha,\quad\mbox{with
}\cos\alpha=\frac{-{(x^{\mathrm{mid}}_{c})}^{2}+{(y^{\mathrm{mid}})}^{2}}{{(x^{\mathrm{mid}}_{c})}^{2}+{(y^{\mathrm{mid}})}^{2}},$
$\displaystyle\mbox{for }c\in\\{\mathrm{lv,rv,sv}\\}.$

It is required that that the total midwall tension at junctions is zero, i.e.,

$f(y^{\mathrm{mid}},V_{\mathrm{sv}}^{\mathrm{mid}}):=\begin{pmatrix}T_{\mathrm{lv}}^{\mathrm{x}}+T_{\mathrm{rv}}^{\mathrm{x}}+T_{\mathrm{sv}}^{\mathrm{x}}\\\
T_{\mathrm{lv}}^{\mathrm{y}}+T_{\mathrm{rv}}^{\mathrm{y}}+T_{\mathrm{sv}}^{\mathrm{y}}\end{pmatrix}\stackrel{{\scriptstyle!}}{{=}}0.$
(43)

Eq. 43 is solved by an iterative Newton scheme

$f^{\prime}(y_{k}^{\mathrm{mid}},V_{k,\mathrm{sv}}^{\mathrm{mid}}){\left(\Delta
y_{k}^{\mathrm{mid}},\Delta
V_{k,\mathrm{sv}}^{\mathrm{mid}}\right)}^{\top}=-f(y_{k}^{\mathrm{mid}},V_{k,\mathrm{sv}}^{\mathrm{mid}}),\quad
k=1,2,\dots$ (44)

and increments $\Delta y_{k}^{\mathrm{mid}}$ and $\Delta
V_{k,\mathrm{sv}}^{\mathrm{mid}}$ are added to $y_{k}^{\mathrm{mid}}$ and
$V_{k,\mathrm{sv}}^{\mathrm{mid}}$. The solution of (44) in the first step,
i.e., for $k=0$ is used to define the ODE updates for the septum

$\dot{V}^{\mathrm{mid}}_{\mathrm{sv}}=\frac{1}{\tau_{\mathrm{sv}}}\Delta
V_{0,\mathrm{sv}}^{\mathrm{mid}},\quad\dot{y}^{\mathrm{mid}}=\frac{1}{\tau_{\mathrm{sv}}}\Delta
y_{0}^{\mathrm{mid}},$ (45)

where $\tau_{\mathrm{sv}}$ is a time constant.

Consequently, the values for the tensions discussed above are updated and the
scheme is iterated until convergence. Midwall volumes are updated by

$\displaystyle V_{c}^{\mathrm{mid}}$
$\displaystyle=V_{c}+\frac{1}{2}\left(V_{c}^{\mathrm{wall}}+V_{\mathrm{sv}}^{\mathrm{wall}}\right),$

long-axis length $l_{c}$ and cross-sectional area $A_{c}$ and of the cavity
are computed by

$\displaystyle l_{c}$
$\displaystyle=2{\left(V_{c}^{\mathrm{mid}}+\frac{1}{2}\left(V_{c}^{\mathrm{wall}}+V_{\mathrm{sv}}^{\mathrm{wall}}\right)\right)}^{1/3},$
(46) $\displaystyle A_{c}$
$\displaystyle=\frac{V_{c}^{\mathrm{mid}}+\frac{1}{20}\left(V_{c}^{\mathrm{wall}}+V_{\mathrm{sv}}^{\mathrm{wall}}\right)}{l_{c}}.$
(47)

Finally, wave impedance $Z_{c}$ is computed according to Eq. 38 and transmural
pressure $p_{c}^{\mathrm{trans}}$ is computed as total axial force

$\displaystyle
p_{c}^{\mathrm{trans}}=2\,\frac{T_{c}^{\mathrm{x}}}{y^{\mathrm{mid}}},\quad\mbox{for
}c\in\\{\mathrm{lv,rv,sv}\\}.$

Assuming the pressure surrounding the ventricular composite to be zero,
internal chamber pressure for the ventricles is now found as

$\displaystyle p_{\mathrm{lv}}$
$\displaystyle=-p_{\mathrm{lv}}^{\mathrm{trans}},$ $\displaystyle
p_{\mathrm{rv}}$ $\displaystyle=p_{\mathrm{rv}}^{\mathrm{trans}}.$

(a)(b)(c)ResistanceValveWall SegmentTubeLV WallSV WallRV
Wall$y$-direction$y^{\mathrm{mid}}$$0$$V_{\mathrm{lv}}^{\mathrm{wall}}$$V_{\mathrm{sv}}^{\mathrm{wall}}$$V_{\mathrm{rv}}^{\mathrm{wall}}$
axis of rotational symmetry intersection juction circle positive $x$-direction
$A_{c}^{\mathrm{mid}}$$\alpha$$V_{c}^{\mathrm{mid}}$$x_{c}^{\mathrm{mid}}$$y^{\mathrm{mid}}$$0$$\frac{1}{C_{c}^{\mathrm{mid}}}$$T_{c}^{\mathrm{mid}}$$T_{c}^{\mathrm{y}}$$T_{c}^{\mathrm{x}}$TriSeg
Model LV SV RV
$V_{\mathrm{lv}}$$p_{\mathrm{lv}}$$V_{\mathrm{rv}}$$p_{\mathrm{rv}}$$V_{\mathrm{rv}}^{\mathrm{wall}}$$V_{\mathrm{sv}}^{\mathrm{wall}}$$V_{\mathrm{lv}}^{\mathrm{wall}}$SystPulm
Right Atrium Left Atrium Syst Veins Pulm Veins Pulm Arteries Syst Arteries
Aortic Valve Pulm Valve Tricuspid Valve Mitral Valve Figure 6: TriSeg model of
ventricular mechanics. (a) The TriSeg model (gray shading) incorporated in the
modular CircAdapt model of the systemic (Syst) and pulmonary (Pulm)
circulations. (b) Cross-section of the ventricular composite. (c) Cross-
section of a single wall segment through the axis of rotational symmetry.
Adapted from [54].

###  A.5 Pericardial mechanics

The four cardiac chambers are supposed to have an additional pressure
component due to the pericardium. Pressure $p_{\mathrm{peri}}$ exerted by the
pericardial sack on atria and ventricles was computed as a non-linear function
of pericardial volume $V_{\mathrm{peri}}$, computed as the sum of blood pool
and wall volumes of the four cardiac chambers:

$\displaystyle V_{\mathrm{peri}}$
$\displaystyle=V_{\mathrm{lv}}+V_{\mathrm{rv}}+V_{\mathrm{la}}+V_{\mathrm{ra}}+V_{\mathrm{lv}}^{\mathrm{wall}}+V_{\mathrm{rv}}^{\mathrm{wall}}+V_{\mathrm{la}}^{\mathrm{wall}}+V_{\mathrm{ra}}^{\mathrm{wall}}$
(48) $\displaystyle p_{\mathrm{peri}}$
$\displaystyle=p_{\mathrm{peri}}^{\mathrm{ref}}{\left(\frac{V_{\mathrm{peri}}}{V_{\mathrm{peri}}^{\mathrm{ref}}}\right)}^{k_{\mathrm{peri}}},$
(49)

where $p_{\mathrm{peri}}^{\mathrm{ref}}$ and
$V_{\mathrm{peri}}^{\mathrm{ref}}$ are constant reference pressure and volume,
respectively, and $k_{\mathrm{peri}}$ defines the degree of non-linearity of
the pressure-volume relation.

Cavity pressures are updated according to

$p_{c}=p_{c}+p_{\mathrm{peri}},\quad\mbox{for
}c\in\\{\mathrm{lv,rv,la,ra}\\}.$ (50)

###  A.6 Periphery

Pulmonary (pulm) and systemic (sys) periphery are modeled as resistances. The
current pressure drop $\Delta p_{py}$, for
$py\in\\{\text{pulm},\text{sys}\\}$, is computed as the difference of the
pressures in the inflow artery $p_{t}^{\mathrm{prox}}$ and the outflow vein
$p_{t}^{\mathrm{dist}}$:

$\Delta p_{py}=p_{t}^{\mathrm{prox}}-p_{t}^{\mathrm{dist}}.$

Using this, the current flow over the periphery is

$q_{py}=q_{py}^{\mathrm{ref}}{\left(r_{py}\frac{\Delta p_{py}}{\Delta
p_{py}^{\mathrm{ref}}}\right)}^{k_{py}},$ (51)

where $\Delta p_{py}^{\mathrm{ref}}$ is the reference arteriovenous pressure
drop; $q_{py}^{\mathrm{ref}}$ is the reference flow over the periphery;
$r_{py}$ is a scaling factor of the ateriovenous resistances; and $k_{py}$ is
a factor that accounts for the nonlinearity of the arteriovenous resistances,
see Tabs. 2 and 3.

###  A.7 Connect modules

Volume change of inflow arteries $\dot{V}_{t}^{\mathrm{prox}}$ and outflow
veins $\dot{V}_{t}^{\mathrm{dist}}$ is now updated by

$\displaystyle\dot{V}_{t}^{\mathrm{dist}}$ $\displaystyle\mathrel{+}=q_{py}$
(52) $\displaystyle\dot{V}_{t}^{\mathrm{prox}}$
$\displaystyle\mathrel{+}=q_{py}$

Computation of time derivative of flow across valves and venous-atrial inlet
requires in input the cross-sectional area of proximal and distal elements to
the channel.

$\displaystyle\dot{V}_{c,t}^{\mathrm{dist}}$ $\displaystyle\mathrel{+}=q_{v}$
(53) $\displaystyle\dot{V}_{c,t}^{\mathrm{prox}}$
$\displaystyle\mathrel{+}=q_{v}$

$\displaystyle p_{c,t}^{\mathrm{prox}}$
$\displaystyle\mathrel{+}=\dot{V}_{c,t}^{\mathrm{prox}}Z_{c,t}^{\mathrm{prox}}$
(54) $\displaystyle p_{c,t}^{\mathrm{dist}}$
$\displaystyle\mathrel{+}=\dot{V}_{c,t}^{\mathrm{dist}}Z_{c,t}^{\mathrm{dist}}$

###  A.8 Valve dynamics

The pressure drop ($\Delta p_{v}$) across a valve is the sum of the effects of
inertia due to acceleration in time and the Bernoulli effect, see [63]

$\Delta
p_{v}=\rho_{\mathrm{b}}\frac{l_{v}}{A_{v}}\dot{q}_{v}+\frac{\rho_{\mathrm{b}}}{2}\left({(v_{v}^{\mathrm{out}})}^{2}-{(v_{v}^{\mathrm{in}})}^{2}\right),$
(55)

where $\rho_{\mathrm{b}}$ is the density of blood, $A_{v}$ is the current
cross-sectional area of the valve, and $l_{v}$ is the length of the channel
with inertia. If not mentioned otherwise this values is estimated as

$l_{v}=\sqrt{A_{v}^{\mathrm{open}}},$

with $A_{v}^{\mathrm{open}}$ the given cross-sectional area of the open valve,
see Tab. 2. For $q_{v}\geq 0$, $v_{v}^{\mathrm{in}}$ is the velocity proximal
to the valve $v_{v}^{\mathrm{prox}}$. $v_{v}^{\mathrm{out}}$ is the maximum of
the blood velocities in the valve region
$v_{v}^{\max}=\max(v_{v}^{\mathrm{dist}},v_{v},v_{v}^{\mathrm{prox}})$. For
$q_{v}<0$ which indicates that the valve is leaking $v_{v}^{\mathrm{in}}$ is
the velocity distal to the valve $v_{v}^{\mathrm{dist}}$ and the outflow
velocity is the maximum of the blood velocities in the valve region
$v_{v}^{\mathrm{out}}=v_{v}^{\max}$. Using $v_{v}=q_{v}/A_{v}$ we can write

$\Delta
p_{v}=p_{v}^{\mathrm{prox}}-p_{v}^{\mathrm{dist}}=\alpha_{v}\dot{q}_{v}+\beta_{v}q_{v}^{2}$
(56)

with

$\alpha_{v}=\rho_{\mathrm{b}}\frac{l_{v}}{A_{v}}$ (57)

the inertia of the channel. The open/closed status of the valve is a function
of pressure drop and flow. Valves are clearly open/closed if both pressure
drop and flow point in the same direction. With forward pressure drop, the
valve opens immediately. With backward pressure and forward flow, the valve is
closing softly by a continuous function

$\displaystyle A_{v}^{\mathrm{closing}}$
$\displaystyle=\sqrt{\frac{x_{v}}{x_{v}^{2}+\Delta
p_{v}^{2}}}\left(A_{v}^{\mathrm{open}}-A_{v}^{\mathrm{leak}}\right)+A_{v}^{\mathrm{leak}}$
(58) $\displaystyle x_{v}$
$\displaystyle=\frac{40\rho_{\mathrm{b}}\,q_{v}\left|q_{v}\right|}{{\left(A_{v}^{\mathrm{open}}\right)}^{2}},$

where $A_{v}^{\mathrm{leak}}$ is the given valve cross-sectional areas of the
closed (regurging) valve. Using this the current cross sectional area of the
valve is

$A_{v}=\begin{cases}A_{v}^{\mathrm{open}}&\ \text{for}\ \Delta p_{v}>0,\\\
A_{v}^{\mathrm{leak}}&\ \text{for}\ \Delta p_{v}<0\ \text{and}\ q_{v}<0,\\\
A_{v}^{\mathrm{closing}}&\ \text{for}\ \Delta p_{v}<0\ \text{and}\ q_{v}>0.\\\
\end{cases}$ (59)

We define

$A_{v}^{\min}=\min\left(A_{v}^{\mathrm{prox}},A_{v},A_{v}^{\mathrm{dist}}\right),$
(60)

with $A_{v}^{\mathrm{prox}}$ and $A_{v}^{\mathrm{dist}}$ the cross-sectional
area of the proximal and distal cavities or tubes respectively, see (20,37,47)
and Fig. 7. Using this $\beta_{v}$ is given as

$\beta_{v}=\begin{cases}\frac{1}{2}\rho_{\mathrm{b}}\left[{\left(\frac{1}{A_{v}^{\min}}\right)}^{2}-{\left(\frac{1}{A_{v}^{\mathrm{prox}}}\right)}^{2}\right]&\
\text{for}\ q_{v}\geq 0,\\\\[10.00002pt]
\frac{1}{2}\rho_{\mathrm{b}}\left[{\left(\frac{1}{A_{v}^{\mathrm{dist}}}\right)}^{2}-{\left(\frac{1}{A_{v}^{\min}}\right)}^{2}\right]&\
\text{for}\ q_{v}<0.\end{cases}$ (61)

Flow over the valve is finally updated using (56) by

$\dot{q}_{v}=\frac{\Delta p_{v}-\beta_{v}q_{v}^{2}}{\alpha_{v}}.$ (62)

(a)(b)$A_{v}^{\mathrm{prox}}$$A_{v}$$A_{v}^{\mathrm{dist}}$$A_{v}^{\mathrm{prox}}$$A_{v}$$A_{v}^{\mathrm{dist}}$
Figure 7: Schematic of the (a) open and (b) regurging valve, based on [64].

###  A.9 Solve ODE system

A Runge–Kutta–Fehlberg method (RKF45), see D, is used to solve the system of
26 ordinary differential equations (ODEs):

  * 8

ODEs: for each of the four tubes and the four cavities we get an ODE to update
the volume using Eqs. 52 and 53.

  * 2

ODEs: for the septum we update midwall volume and the radius according to
(45).

  * 10

ODEs: for the sarcomeres of each cavity and the septum we update sarcomere
contracting length and contractility using (26–27).

  * 6

ODEs: for each of the four valves and the two outlets we update flow by (62).

Parameter |  Value | Unit |  Description
---|---|---|---
_General_
$\rho_{\mathrm{b}}$ |  $1050.0$ | $\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{{}^{3}}$ |  blood density
$t_{\mathrm{cycle}}$ |  $0.585$ | $\mathrm{s}$ |  cycle time ($=1/\mathrm{heartrate}$)
_Tubes:_ aorta (ao), arteria pulmonalis (ap), venae cavae (vc), and venae
pulmonales (vp)
$A_{\mathrm{t}}^{\mathrm{wall}}$ |  274 (ao), 141 (ap), 58 (vc), 85 (vp) | $\mathrm{m}\mathrm{m}^{2}$ |  cross-sectional wall area
$l_{\mathrm{t}}$ |  500 (ao), 400 (vc), 200 (ap, vp) | $\mathrm{mm}$ |  length of vessel
$A_{\mathrm{t}}^{\mathrm{ref}}$ |  adjacent valve area | $\mathrm{m}\mathrm{m}^{2}$ |  initial cross sectional area
$k_{t}$ |  5 (ao), 8 (ap), 10 (vc, vp) | [-] |  stiffness exponent
_Chambers:_ left (lv) and right (rv) ventricle; left (la) and right atrium
(ra)
$V_{\mathrm{c}}$ |  57.0 (lv), 75.3 (rv), 44.2 (la), 54.4 (ra) | $\mathrm{mL}$ |  cavity volume
$h_{\mathrm{c}}^{\mathrm{wall}}$ |  15.0 (lv), 4.0 (rv), 2.0 (la), 2.0 (ra) | $\mathrm{mm}$ |  constant wall thickness
$\Delta t_{\mathrm{c}}^{\mathrm{act}}$ |  0.1 (lv), 0.1 (rv), 0.02 (la), 0.0 (ra) | $\mathrm{s}$ |  delays of onset of activation in each beat starting at $t^{\mathrm{bt}}$; $t_{c}^{\mathrm{act}}=t^{\mathrm{bt}}+\Delta t_{\mathrm{c}}^{\mathrm{act}}$
_Valves:_ aortic (av), pulmonary (pv), mitral (mv), and tricuspid (tv) valve;
 _Valves:_ pulmonary (po) and systemic (so) outlet
$A_{v}^{\mathrm{open}}$ |  500 (mv, tv) 400 (av, pv, so, po) | $\mathrm{m}\mathrm{m}^{2}$ |  valve cross-sectional area
$A_{v}^{\mathrm{leak}}$ |  0 (av, pv, mv, tv) $A_{v}^{\mathrm{open}}$ (so, po) | $\mathrm{m}\mathrm{m}^{2}$ |  cross-sectional area of closed/regurging valve
_Periphery:_ systemic (sys) and pulmonary (pulm) circulation
$\Delta p_{py}^{\mathrm{ref}}$ |  1.5 (pulm) 10.0 (sys) | $\mathrm{kPa}$ |  blood pressure drop in pulmonary/systemic circulation
$q_{py}^{\mathrm{ref}}$ |  $85$ (pulm, sys) | $\mathrm{m}\mathrm{L}\mathrm{/}\mathrm{s}$ |  reference pulmonary/systemic flow
$r_{py}$ |  1 (pulm) 2 (sys) | [-] |  resistance scaling factor
Table 2: Input parameters for the CircAdapt model. Adjusted to match patient-specific data. Parameter |  Value | Unit |  Description
---|---|---|---
_Tubes_
$p_{t}^{\mathrm{ref}}$ |  12.0 (ao) 0.5 (vp) 0.12 (vc) 1.8 (ap) | $\mathrm{kPa}$ |  reference tube pressure
_Sarcomeres_
$L^{\mathrm{s,ref}}$ |  2.00 | $\mathrm{\SIUnitSymbolMicro m}$ |  reference sarcomere length
$L^{\mathrm{elast,iso}}$ |  0.04 | $\mathrm{\SIUnitSymbolMicro m}$ |  length of isometrically stressed series elastic element
$v^{\max}$ |  7 (lv, rv) 14 (la, ra) | $\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{/}\mathrm{s}$ |  reference shortening velocity
$\tau^{\mathrm{R}}$ |  0.25 (lv, rv) 0.4 (la, ra) | [-] |  ratio rise time to $t_{c}^{\mathrm{act,ref}}$
$\tau^{\mathrm{D}}$ |  0.25 (lv, rv) 0.4 (la, ra) | [-] |  ratio decay time to $t_{c}^{\mathrm{act,ref}}$
$L^{\mathrm{act0,ref}}$ |  1.51 | $\mathrm{\SIUnitSymbolMicro m}$ |  contractile element length with zero active stress
$L^{\mathrm{pas0,ref}}$ |  1.80 | $\mathrm{\SIUnitSymbolMicro m}$ |  sarcomere length with zero passive stress
$\sigma^{\mathrm{act,\max}}$ |  120 (lv, rv) 84 (la, ra) | $\mathrm{kPa}$ |  maximal isometric stress
$\sigma^{\mathrm{pas,\max}}$ |  22 (lv, rv) 50 (la, ra) | $\mathrm{kPa}$ |  maximal passive stress
$\,\mathrm{d}L^{\mathrm{s,pas}}$ |  0.6 | $\mathrm{\SIUnitSymbolMicro m}$ |
_TriSeg Module_
$\tau_{\mathrm{sv}}$ |  0.005 | [-] |  time constant
_Pericardium_
$p_{\mathrm{peri}}^{\mathrm{ref}}$ |  0.005 | [-] |  constant reference pressure
$V_{\mathrm{peri}}^{\mathrm{ref}}$ |  0.005 | [-] |  constant reference volume
_Periphery_
$k_{py}$ |  2 (pulm) 1 (sys) | [-] |  nonlinearity exponent
Table 3: Default parameters for the CircAdapt model, fitted to general
experimental data in [38].

##  Appendix B Finite Element Formulation

###  B.1 Variational Formulation

We first ignore the acceleration term in Eq. 3 and look at the stationary
version of the boundary value problem (3–4) and (15). For the full nonlinear
elastodynamics problem see C. The stationary boundary value problem is
formally equivalent to the equations

$\displaystyle\langle\mathcal{A}_{0}({\mathbf{u}}),{\mathbf{v}}\rangle_{\Omega_{0}}-\langle\mathcal{F}_{0}({\mathbf{u}},p_{c}),{\mathbf{v}}\rangle_{\Omega_{0}}$
$\displaystyle={\mathbf{0}},$ (63) $\displaystyle\langle
V_{c}^{\mathrm{PDE}}({\mathbf{u}}),q\rangle_{\Omega_{0}}-\langle
V_{c}^{\mathrm{ODE}}(p_{c}),q\rangle_{\Omega_{0}}$ $\displaystyle=0,$ (64)

which is valid for all smooth enough vector fields ${\mathbf{v}}$ vanishing on
the Dirichlet boundary ${\Gamma}_{0,\mathrm{D}}$, testfunctions $q$ that are
$1$ for the cavity $c$ and $0$ otherwise, the duality pairing
$\langle\cdot,\cdot\rangle_{\Omega_{0}}$, and cavities
$c\in\\{\mathrm{lv,rv,la,ra}\\}$. The second term on the left hand side of the
variational equation (63) has the physical interpretation of the rate of
internal mechanical work and is given by

$\langle\mathcal{A}_{0}({\mathbf{u}}),{\mathbf{v}}\rangle_{\Omega_{0}}:=\int_{\Omega_{0}}\bm{\mathbf{S}}({\mathbf{u}}):\bm{\mathbf{\Sigma}}({\mathbf{u}},{\mathbf{v}})\,\mathrm{d}{\mathbf{X}},$
(65)

with the second Piola–Kirchhoff stress tensor $\bm{\mathbf{S}}$, see (5), and
the directional derivative of the Green–Lagrange strain tensor
$\bm{\mathbf{\Sigma}}({\mathbf{u}},{\mathbf{v}})$, see [35, 65]. The weak form
of the contribution of pressure loads (63), right term, is computed using (4)

$\langle\mathcal{F}_{0}({\mathbf{u}},p_{c}),{\mathbf{v}}\rangle_{\Omega_{0}}=-p_{c}\int\limits_{\Gamma_{0,\mathrm{N}}}J\,\bm{\mathbf{F}}^{-\top}({\mathbf{u}})\,{\mathbf{n}^{\mathrm{out}}_{0}}\cdot{\mathbf{v}}\,\,\mathrm{d}s_{\bm{\mathbf{X}}}.$
(66)

The first term of the coupling equation (64) is computed from (14) using
Nanson’s formula and ${\mathbf{x}}=\bm{\mathbf{X}}+{\mathbf{u}}$ by

$\langle
V_{c}^{\mathrm{PDE}}({\mathbf{u}}),q\rangle_{\Omega_{0}}=\frac{1}{3}\int\limits_{\Gamma_{0,\mathrm{N}}}(\bm{\mathbf{X}}+{\mathbf{u}})\cdot
J\bm{\mathbf{F}}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\,q\,\mathrm{d}s_{\bm{\mathbf{X}}},$
(67)

The second term of (64) is computed using the lumped CircAdapt model, see 2.3,
for $c\in\\{\mathrm{lv,rv,la,ra}\\}$.

###  B.2 Consistent Linearization

To solve the nonlinear variational equations (63)–(64), with a FE approach we
first apply a Newton–Raphson scheme, see [66]. Given a nonlinear and
continuously differentiable operator $F\colon X\to Y$ a solution to $F(x)=0$
can be approximated by

$\displaystyle x^{k+1}$ $\displaystyle=x^{k}+\Delta x,$
$\displaystyle\left.\frac{\partial F}{\partial x}\right|_{x=x^{k}}\Delta x$
$\displaystyle=-F(x^{k}),$

which is looped until convergence. In our case, we have
$X=V_{{\mathbf{g}}_{D}}\times\mathbb{R}$, $Y=\mathbb{R}^{2}$, $\Delta
x={(\Delta{\mathbf{u}},\Delta p_{c})}^{\top}$,
$x^{k}={({\mathbf{u}}^{k},p_{c}^{k})}^{\top}$, and
$F={(R_{\mathbf{u}},R_{\mathrm{p}})}^{\top}$. We obtain the following
linearized saddle-point problem for each $({\mathbf{u}}^{k},p_{c}^{k})\in
V_{{\mathbf{g}}_{D}}\times\mathbb{R}$, find $(\Delta{\mathbf{u}},\Delta
p_{c})\in V_{{\mathbf{0}}}\times\mathbb{R}$ such that

$\displaystyle\langle\Delta{\mathbf{u}},A_{0}^{\prime}({\mathbf{u}}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}+\langle\Delta{\mathbf{u}},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}$
$\displaystyle+\langle\Delta
p_{c},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}$
$\displaystyle=-\langle
R_{\mathbf{u}}({\mathbf{u}}^{k},p_{c}^{k}),{\mathbf{v}}\rangle_{\Omega_{0}},$
(68)
$\displaystyle\langle\Delta{\mathbf{u}},V_{c}^{\mathrm{PDE}}({\mathbf{u}}^{k})\,q\rangle_{\Omega_{0}}-\langle\Delta
p_{c},V_{c}^{\mathrm{ODE}}(p_{c}^{k})\,q\rangle_{\Omega_{0}}$
$\displaystyle=-\langle
R_{\mathrm{p}}({\mathbf{u}}^{k},p_{c}^{k}),q\rangle_{\Omega_{0}},$ (69)

with the updates

$\displaystyle{\mathbf{u}}^{k+1}$
$\displaystyle={\mathbf{u}}^{k}+\Delta{\mathbf{u}},$ (70) $\displaystyle
p_{c}^{k+1}$ $\displaystyle=p_{c}^{k}+\Delta p_{c},$ (71)

and the particular terms are introduced below. The Gâteaux derivative of (63)
with respect to the displacement change update $\Delta{\mathbf{u}}$ yields the
first

$\displaystyle\langle\Delta{\mathbf{u}},A_{0}^{\prime}({\mathbf{u}}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}:$
$\displaystyle=\left.D_{\Delta{\mathbf{u}}}\langle\mathcal{A}_{0}({\mathbf{u}}),{\mathbf{v}}\rangle_{\Omega_{0}}\right|_{{\mathbf{u}}={\mathbf{u}}^{k}}$
$\displaystyle=\int\limits_{\Omega_{0}}\bm{\mathbf{S}}_{k}:\bm{\mathbf{\Sigma}}(\Delta{\mathbf{u}},{\mathbf{v}})\,\mathrm{d}\bm{\mathbf{X}}+\int\limits_{\Omega_{0}}\bm{\mathbf{\Sigma}}({\mathbf{u}}^{k},\Delta{\mathbf{u}}):\mathbb{C}_{k}:\bm{\mathbf{\Sigma}}({\mathbf{u}}^{k},{\mathbf{v}})\,\mathrm{d}\bm{\mathbf{X}},$
(72)

and second term of (B.2)

$\displaystyle\langle\Delta{\mathbf{u}},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}:$
$\displaystyle=\left.D_{\Delta{\mathbf{u}}}\langle\mathcal{F}_{0}({\mathbf{u}},p_{c}),{\mathbf{v}}\rangle_{\Omega_{0}}\right|_{{\mathbf{u}}={\mathbf{u}}^{k},p_{c}=p_{c}^{k}}$
$\displaystyle=p^{k}_{c}\int_{\Gamma_{0,\mathrm{N}}}J_{k}\bm{\mathbf{F}}_{k}^{-\top}\operatorname{Grad}^{\top}\\!\Delta{\mathbf{u}}\,\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\cdot{\mathbf{v}}\,\,\mathrm{d}s_{\bm{\mathbf{X}}}$
$\displaystyle-p^{k}_{c}\int_{\Gamma_{0,\mathrm{N}}}J_{k}(\bm{\mathbf{F}}_{k}^{-\top}:\operatorname{Grad}\Delta{\mathbf{u}})\,\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\cdot{\mathbf{v}}\,\,\mathrm{d}s_{\bm{\mathbf{X}}},$
(73)

with abbreviations

$\displaystyle\bm{\mathbf{F}}_{k}$
$\displaystyle:=\bm{\mathbf{F}}({\mathbf{u}}^{k}),\
J_{k}:=\det(\bm{\mathbf{F}}^{k}),\
\bm{\mathbf{S}}_{k}:=\left.\bm{\mathbf{S}}\right|_{{\mathbf{u}}={\mathbf{u}}^{k}},\
\mathbb{C}_{k}:=\left.\mathbb{C}\right|_{{\mathbf{u}}={\mathbf{u}}^{k}}.$

The Gâteaux derivative of (63) with respect to the pressure change update
$\Delta p_{c}$ yields the third term of (B.2)

$\displaystyle\langle\Delta
p_{c},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\,{\mathbf{v}}\rangle_{\Omega_{0}}:$
$\displaystyle=\left.D_{\Delta
p_{c}}\langle\mathcal{F}_{0}({\mathbf{u}},p_{c}),{\mathbf{v}}\rangle_{\Omega_{0}}\right|_{{\mathbf{u}}={\mathbf{u}}^{k},p_{c}=p_{c}^{k}}$
$\displaystyle=-\Delta
p_{c}\int_{\Gamma_{0,\mathrm{N}}}J_{k}\,\bm{\mathbf{F}}_{k}^{-\top}\,{\mathbf{n}^{\mathrm{out}}_{0}}\cdot{\mathbf{v}}\,\,\mathrm{d}s_{{\mathbf{X}}}.$
(74)

The residual $R_{\mathbf{u}}$, i.e., the right hand side of (B.2), is computed
as

$\langle
R_{\mathbf{u}}({\mathbf{u}}^{k},p_{c}^{k}),{\mathbf{v}}\rangle_{\Omega_{0}}:=\langle
A_{0}({\mathbf{u}}^{k}),{\mathbf{v}}\rangle_{\Omega_{0}}-\langle{\mathbf{\mathcal{F}}}_{0}({\mathbf{u}}^{k},p_{c}^{k}),{\mathbf{v}}\rangle_{\Omega_{0}}.$
(75)

From (67), using the known relations, see, e.g., [65],

$\displaystyle\frac{\partial
J}{\partial\bm{\mathbf{F}}}:\operatorname{Grad}\Delta{\mathbf{u}}$
$\displaystyle=J\bm{\mathbf{F}}^{-\top}:\operatorname{Grad}\Delta{\mathbf{u}}$
$\displaystyle\frac{\partial\bm{\mathbf{F}}^{-\top}}{\partial\bm{\mathbf{F}}}:\operatorname{Grad}\Delta{\mathbf{u}}$
$\displaystyle=-\bm{\mathbf{F}}^{-\top}{(\operatorname{Grad}\Delta{\mathbf{u}})}^{\top}\bm{\mathbf{F}}^{-\top}$

we can calculate the first term of (69) as the Gâteaux derivative with respect
to the update $\Delta{\mathbf{u}}$

$\displaystyle\langle\Delta{\mathbf{u}},V_{c}^{\mathrm{PDE}}({\mathbf{u}}^{k})\,q\rangle_{\Omega_{0}}:$
$\displaystyle=\left.D_{\Delta{\mathbf{u}}}\langle
V_{c}^{\mathrm{PDE}}({\mathbf{u}}),q\rangle_{\Omega_{0}}\right|_{{\mathbf{u}}={\mathbf{u}}^{k}}$
$\displaystyle=D_{\Delta{\mathbf{u}}}\frac{1}{3}\int_{\Gamma_{0,\mathrm{N}}}\left({\mathbf{X}}+{\mathbf{u}}^{k}\right)\cdot
J_{k}\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\,q\,\mathrm{d}s_{\bm{\mathbf{X}}}$
$\displaystyle=\frac{1}{3}\int_{\Gamma_{0,\mathrm{N}}}J_{k}(\bm{\mathbf{F}}_{k}^{-\top}:\operatorname{Grad}\Delta{\mathbf{u}}){\mathbf{x}}\cdot\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\,q\,\mathrm{d}s_{\bm{\mathbf{X}}}$
$\displaystyle\quad-\frac{1}{3}\int_{\Gamma_{0,\mathrm{N}}}J_{k}{\mathbf{x}}\cdot\bm{\mathbf{F}}_{k}^{-\top}{(\operatorname{Grad}\Delta{\mathbf{u}})}^{\top}\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\,q\,\mathrm{d}s_{\bm{\mathbf{X}}}$
$\displaystyle\quad+\frac{1}{3}\int_{\Gamma_{0,\mathrm{N}}}J_{k}\Delta{\mathbf{u}}\cdot\bm{\mathbf{F}}_{k}^{-\top}{\mathbf{n}^{\mathrm{out}}_{0}}\,q\,\mathrm{d}s_{\bm{\mathbf{X}}},$
(76)

with $q$ a testfunction that is $1$ for the surface of cavity $c$,
$\Gamma_{0,c}$, and $0$ otherwise.

The second term of (69) is computed as a numerical derivative

$\displaystyle\langle\Delta
p_{c},V_{c}^{\mathrm{ODE}}(p_{c}^{k})\,q\rangle_{\Omega_{0}}:$
$\displaystyle=\left.D_{\Delta p_{c}}\langle
V_{c}^{\mathrm{ODE}}(p_{c})\,q\rangle_{\Omega_{0}}\right|_{p_{c}=p_{c}^{k}}$
$\displaystyle=\frac{1}{\epsilon}\left(V_{c}^{\mathrm{ODE}}(p_{c}^{k}+\epsilon)-V_{c}^{\mathrm{ODE}}(p_{c}^{k})\right)q,$
(77)

where $\epsilon=p_{c}^{k}\sqrt{\epsilon_{\mathrm{m}}}$ is chosen according to
[67, Chapter 5.7] with $\epsilon_{\mathrm{m}}=2^{-52}\approx 2.2*10^{-16}$ the
machine accuracy.

Finally, the residual $R_{\mathrm{p}}$, i.e., the right hand side of (69), is
computed as

$\langle
R_{\mathrm{p}}({\mathbf{u}}^{k},p_{c}^{k}),q\rangle_{\Omega_{0}}:=\langle
V_{c}^{\mathrm{PDE}}({\mathbf{u}}),q\rangle_{\Omega_{0}}-\langle
V_{c}^{\mathrm{ODE}}(p_{c}),q\rangle_{\Omega_{0}}.$ (78)

###  B.3 Assembling of the block matrices

To apply the finite element method (FEM) we consider an admissible
decomposition of the computational domain $\Omega\subset\mathbb{R}^{3}$ into
$M$ tetrahedral elements $\tau_{j}$ and introduce a conformal finite element
space

$\displaystyle X_{h}\subset H^{1}(\Omega_{0}),\ N=\text{dim}\,X_{h}$

of piecewise polynomial continuous basis functions $\varphi_{i}$. The
linearized variational problem (B.2)–(69) and a Galerkin FE discretization
result in solving the block system to find
$\delta\underline{u}\in\mathbb{R}^{3N}$ and
$\delta\underline{p}_{\mathtt{c}}\in\mathbb{R}^{N_{\mathrm{cav}}}$ such that

$\bm{\mathbf{K}}^{\prime}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\begin{pmatrix}\delta\underline{u}\\\
\delta\underline{p}_{\mathtt{c}}\end{pmatrix}=-\underline{K}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k}),\quad\underline{K}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k}):=-\begin{pmatrix}\underline{R}_{{\mathbf{u}}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\\\
\underline{R}_{p}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\end{pmatrix},$

i.e.,

$\displaystyle\begin{pmatrix}(\bm{\mathbf{A}}^{\prime}-\bm{\mathbf{M}}^{\prime})(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})&\bm{\mathbf{B}}^{\prime}_{\mathrm{p}}(\underline{u}^{k})\\\
\bm{\mathbf{B}}^{\prime}_{\mathbf{u}}(\underline{u}^{k})&\bm{\mathbf{C}}^{\prime}(\underline{p}_{\mathtt{c}}^{k})\\\
\end{pmatrix}\begin{pmatrix}\delta\underline{u}\\\
\delta\underline{p}_{\mathtt{c}}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}\underline{A}(\underline{u}^{k})-\underline{B}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\\\
\underline{V}_{c}^{\mathrm{PDE}}(\underline{u}^{k})-\underline{V}_{c}^{\mathrm{ODE}}(\underline{p}_{\mathtt{c}}^{k})\end{pmatrix},$
(79) $\displaystyle\underline{u}^{k+1}$
$\displaystyle=\underline{u}^{k}+\delta\underline{u},$ (80)
$\displaystyle\underline{p}_{\mathtt{c}}^{k+1}$
$\displaystyle=\underline{p}_{\mathtt{c}}^{k}+\delta\underline{p}_{\mathtt{c}}$
(81)

with the solution vectors $\underline{u}^{k}\in\mathbb{R}^{3N}$ and
$\underline{p}_{\mathtt{c}}^{k}\in\mathbb{R}^{N_{\mathrm{cav}}}$ at the $k$-th
Newton step. The tangent stiffness matrix
$\bm{\mathbf{A}}^{\prime}\in\mathbb{R}^{3N\times 3N}$ is calculated from (72)
according to

$\bm{\mathbf{A}}^{\prime}(\underline{u}^{k})[j,i]:=\langle\bm{\varphi}_{i},\mathcal{A}_{0}^{\prime}({\mathbf{u}}^{k})\,\bm{\varphi}_{j}\rangle_{\Omega_{0}}$
(82)

and the mass matrix $\bm{\mathbf{M}}^{\prime}\in\mathbb{R}^{3N\times 3N}$ is
calculated from (73) according to

$\bm{\mathbf{M}}^{\prime}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})[j,i]:=\langle\bm{\varphi}_{i},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\,\bm{\varphi}_{j}\rangle_{\Omega_{0}},$
(83)

see also [35, 65].

The off-diagonal matrices
$\bm{\mathbf{B}}^{\prime}_{{\mathbf{u}}}\in\mathbb{R}^{3N\times
N_{\mathrm{cav}}}$ and
$\bm{\mathbf{B}}^{\prime}_{\mathrm{p}}\in\mathbb{R}^{N_{\mathrm{cav}}\times
3N}$ in (79) are assembled using (76)

$\bm{\mathbf{B}}^{\prime}_{{\mathbf{u}}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})[i,j]=\langle\bm{\varphi}_{j},V_{c}^{\mathrm{PDE}}({\mathbf{u}}^{k})\hat{\varphi_{i}}\rangle_{\Omega_{0}},\quad
i=1,\dots,N_{\mathrm{cav}}$ (84)

and using (74)

$\bm{\mathbf{B}}^{\prime}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})[i,j]=\langle\hat{\varphi_{j}},\mathcal{F}_{0}^{\prime}({\mathbf{u}}^{k},p_{c}^{k})\bm{\varphi}_{i}\rangle_{\Omega_{0}},\quad
j=1,\dots,N_{\mathrm{cav}},$ (85)

with the constant shape function $\hat{\varphi_{j}}=1$ if
$\tau_{j}\in\Gamma_{0,c}$ and $\hat{\varphi_{j}}=0$ if
$\tau_{j}\notin\Gamma_{0,c}$ for $c\in\\{\mathrm{lv,rv,la,ra}\\}$.

Using a technique as described in [56, Sect. 4.2] this assembling procedure
can be simplified for closed cavities such that

$\bm{\mathbf{B}}^{\prime}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})={\left[\bm{\mathbf{B}}^{\prime}_{{\mathbf{u}}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})\right]}^{\top}.$

The circulatory compliance matrix
$\bm{\mathbf{C}}^{\prime}(\underline{p}_{\mathtt{c}}^{k})\in\mathbb{R}^{N_{\mathrm{cav}}\times
N_{\mathrm{cav}}}$ is computed from (77) as

$\bm{\mathbf{C}}^{\prime}(\underline{p}_{\mathtt{c}}^{k})[i,j]=\langle\hat{\varphi_{j}},V_{c}^{\mathrm{ODE}}(p_{c}^{k})\,\hat{\varphi_{i}}\rangle_{\Omega_{0}},\quad
i,j=1,\dots,N_{\mathrm{cav}},$ (86)

with the constant shape function $\hat{\varphi_{i}},\hat{\varphi_{j}}=1$ for
cavity c and 0 otherwise, leading to a diagonal matrix.

The terms on the upper right hand side $\underline{A}\in\mathbb{R}^{3N}$,
$\underline{B}_{\mathrm{p}}\in\mathbb{R}^{3N}$ are constructed using (75)
resulting in
$\underline{R}_{{\mathbf{u}}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})=\underline{A}(\underline{u}^{k})-\underline{B}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})$
with

$\underline{A}(\underline{u}^{k})[i]:=\langle\mathcal{A}_{0}({\mathbf{u}}^{k}),\bm{\varphi}_{i}\rangle_{\Omega_{0}}$
(87)

and

$\underline{B}_{\mathrm{p}}(\underline{u}^{k},\underline{p}_{\mathtt{c}}^{k})[i]:=\langle\mathcal{F}_{0}({\mathbf{u}}^{k},p_{c}^{k}),\bm{\varphi}_{i}\rangle_{\Omega_{0}}.$
(88)

Finally, the lower right hand side in (79),
$\underline{R}_{p}(\underline{u}^{k},\underline{p}_{\mathrm{c}}^{k})=\underline{V}^{\mathrm{PDE}}(\underline{u}^{k})-\underline{V}^{\mathrm{ODE}}(\underline{p}_{\mathrm{c}}^{k})\in\mathbb{R}^{N_{\mathrm{cav}}}$,
is assembled from (78) with

$\underline{V}^{\mathrm{PDE}}(\underline{u}^{k})[i]=\langle
V_{c}^{\mathrm{PDE}}({\mathbf{u}}),\hat{\varphi_{i}}\rangle_{\Omega_{0}},\quad
i=1,\dots,N_{\mathrm{cav}},$ (89)

and

$\underline{V}^{\mathrm{ODE}}(\underline{p}_{\mathrm{c}}^{k})[i]=\langle
V_{c}^{\mathrm{ODE}}(p_{c}),\hat{\varphi_{i}}\rangle_{\Omega_{0}},\quad
i=1,\dots,N_{\mathrm{cav}}.$ (90)

##  Appendix C Generalized-$\alpha$ Scheme

After standard discretization we rewrite Eq. 3 using Eqs. 87 and 88 as a
nonlinear ODE reading

$\rho_{0}\bm{\mathbf{M}}_{\alpha}\ddot{\underline{u}}(t)+\underline{R}_{{\mathbf{u}}}(\underline{u},t)=\underline{0},$
(91)

with the mass matrix

$\displaystyle\bm{\mathbf{M}}_{\alpha}[i,j]$
$\displaystyle:=\int_{\Omega_{0}}\bm{\varphi}_{i}({\mathbf{X}})\cdot\bm{\varphi}_{j}({\mathbf{X}})\,\mathrm{d}{\mathbf{X}}.$

Following [68] we reformulate Eq. 91 as a first order ODE system by
introducing the velocity $\underline{v}$

$\displaystyle\rho_{0}\bm{\mathbf{M}}_{\alpha}\dot{\underline{v}}(t)$
$\displaystyle+\underline{R}_{{\mathbf{u}}}(\underline{u},t)$
$\displaystyle=\underline{0},$ (92)
$\displaystyle\bm{\mathbf{M}}_{\alpha}\dot{\underline{u}}(t)$
$\displaystyle-\bm{\mathbf{M}}_{\alpha}\underline{v}(t)$
$\displaystyle=\underline{0}$ (93)

and apply a generalized-$\alpha$ approach [69]. To this end we define three
parameters

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

where the _spectral radius_ $\rho_{\infty}$ is a parameter between $0$ and
$1$. With this we introduce

$\displaystyle\dot{\underline{v}}_{n+\alpha_{\mathrm{m}}}$
$\displaystyle:=\alpha_{\mathrm{m}}$ $\displaystyle\dot{\underline{v}}_{n+1}$
$\displaystyle+(1-\alpha_{\mathrm{m}})$
$\displaystyle\dot{\underline{v}}_{n},$
$\displaystyle\dot{\underline{u}}_{n+\alpha_{\mathrm{m}}}$
$\displaystyle:=\alpha_{\mathrm{m}}$ $\displaystyle\dot{\underline{u}}_{n+1}$
$\displaystyle+(1-\alpha_{\mathrm{m}})$
$\displaystyle\dot{\underline{u}}_{n},$
$\displaystyle\underline{v}_{n+\alpha_{\mathrm{f}}}$
$\displaystyle:=\alpha_{\mathrm{f}}$ $\displaystyle\underline{v}_{n+1}$
$\displaystyle+(1-\alpha_{\mathrm{f}})$ $\displaystyle\underline{v}_{n},$
$\displaystyle\underline{u}_{n+\alpha_{\mathrm{f}}}$
$\displaystyle:=\alpha_{\mathrm{f}}$ $\displaystyle\underline{u}_{n+1}$
$\displaystyle+(1-\alpha_{\mathrm{f}})$ $\displaystyle\underline{u}_{n},$

and reformulate Eq. 92 as

$\displaystyle\rho_{0}\bm{\mathbf{M}}_{\alpha}\dot{\underline{v}}_{n+\alpha_{\mathrm{m}}}$
$\displaystyle+\underline{R}_{{\mathbf{u}}}(\underline{u}_{n+\alpha_{\mathrm{f}}})$
$\displaystyle=\underline{0},$ (94)
$\displaystyle\bm{\mathbf{M}}_{\alpha}\dot{\underline{u}}_{n+\alpha_{\mathrm{m}}}$
$\displaystyle-\bm{\mathbf{M}}_{\alpha}\underline{v}_{n+\alpha_{\mathrm{f}}}$
$\displaystyle=\underline{0}.$ (95)

Here, the second equation gives us

$\displaystyle\dot{\underline{u}}_{n+\alpha_{\mathrm{m}}}=\underline{v}_{n+\alpha_{\mathrm{f}}}$

and we get for the velocity update

$\displaystyle\underline{v}_{n+1}=\frac{\alpha_{\mathrm{m}}}{\alpha_{\mathrm{f}}\gamma\Delta
t}\left(\underline{u}_{n+1}-\underline{u}_{n}\right)+\frac{\gamma-\alpha_{\mathrm{m}}}{\gamma\alpha_{\mathrm{f}}}\dot{\underline{u}}_{n}+\frac{\alpha_{\mathrm{f}}-1}{\alpha_{\mathrm{f}}}\underline{v}_{n}.$
(96)

From this and the relationship by Newmark [70]

$\displaystyle\underline{u}_{n+1}$ $\displaystyle=\underline{u}_{n}$
$\displaystyle+\Delta
t\left(\gamma\dot{\underline{u}}_{n+1}+(1-\gamma)\dot{\underline{u}}_{n}\right),$
$\displaystyle\underline{v}_{n+1}$ $\displaystyle=\underline{v}_{n}$
$\displaystyle+\Delta
t\left(\gamma\dot{\underline{v}}_{n+1}+(1-\gamma)\dot{\underline{v}}_{n}\right),$

we obtain

$\displaystyle\dot{\underline{v}}_{n+1}=\frac{\alpha_{\mathrm{m}}}{\alpha_{\mathrm{f}}\gamma^{2}\Delta
t^{2}}\left(\underline{u}_{n+1}-\underline{u}_{n}\right)-\frac{1}{\alpha_{\mathrm{f}}\gamma\Delta
t}\underline{v}_{n}+\frac{\gamma-1}{\gamma}\dot{\underline{v}}_{n}+\frac{\gamma-\alpha_{\mathrm{m}}}{\alpha_{\mathrm{f}}\gamma^{2}\Delta
t}.$ (97)

Hence, we can rewrite the whole first order system only dependent on the
unknowns $\underline{u}_{n+1}$.

#### Newton’s method for the generalized-$\alpha$ scheme

For the implementation of Newton’s method we compute

$\frac{\partial\dot{\underline{v}}_{n+\alpha_{\mathrm{m}}}}{\partial\underline{u}_{n+1}}=\frac{\alpha_{\mathrm{m}}^{2}}{\alpha_{\mathrm{f}}\gamma^{2}\Delta
t^{2}},\quad\frac{\partial\underline{v}_{n+\alpha_{\mathrm{f}}}}{\partial\underline{u}_{n+1}}=\frac{\alpha_{\mathrm{f}}\alpha_{\mathrm{m}}}{\alpha_{\mathrm{f}}\gamma\Delta
t}=\frac{\alpha_{\mathrm{m}}}{\gamma\Delta
t},\quad\frac{\partial\underline{u}_{n+\alpha_{\mathrm{f}}}}{\partial\underline{u}_{n+1}}=\alpha_{\mathrm{f}}.$
(98)

To calculate the solution at the current timestep we assume that we know
$\underline{u}_{n}$, $\dot{\underline{u}}_{n}$, $\underline{v}_{n}$ and
$\dot{\underline{v}}_{n}$ from the previous time step $n$ and get from Eq. 94
for the residual

$\underline{R}_{\alpha}(\underline{u}_{n+1}^{k}):=-\rho_{0}\bm{\mathbf{M}}_{\alpha}\dot{\underline{v}}^{k}_{n+\alpha_{\mathrm{m}}}-\underline{R}_{\mathbf{u}}(\underline{u}^{k}_{n+\alpha_{\mathrm{f}}}),$
(99)

with
$\dot{\underline{v}}^{k}_{n+\alpha_{\mathrm{m}}}:=\dot{\underline{v}}_{n+\alpha_{\mathrm{m}}}(\underline{u}^{k}_{n+1})$
and
$\underline{u}^{k}_{n+\alpha_{\mathrm{f}}}:=\underline{u}_{n+\alpha_{\mathrm{f}}}(\underline{u}^{k}_{n+1})$.
To avoid numerical dissipation we add the following mass and stiffness damping
terms

$\displaystyle\bm{\mathbf{D}}_{\mathrm{mass}}(\underline{u}_{n+1}^{k})$
$\displaystyle=\rho_{0}$ $\displaystyle\beta_{\mathrm{mass}}$
$\displaystyle\bm{\mathbf{M}}_{\alpha}$
$\displaystyle\underline{v}^{k}_{n+\alpha_{\mathrm{f}}},$ (100)
$\displaystyle\bm{\mathbf{D}}_{\mathrm{stiff}}(\underline{u}_{n+1}^{k})$
$\displaystyle=$ $\displaystyle\beta_{\mathrm{stiff}}$
$\displaystyle\bm{\mathbf{A}}^{\prime}$
$\displaystyle\underline{v}^{k}_{n+\alpha_{\mathrm{f}}},$ (101)

to the residual (99) with
$\underline{v}^{k}_{n+\alpha_{\mathrm{f}}}:=\underline{v}_{n+\alpha_{\mathrm{f}}}(\underline{u}^{k}_{n+1})$.

The tangent stiffness matrix is now calculated using (98) as

$\displaystyle\bm{\mathbf{A}}^{\prime}_{\alpha}(\underline{u}^{k}_{n+1},\underline{p}_{\mathtt{c},{n+1}}^{k})$
$\displaystyle:=\rho_{0}\frac{\partial\dot{\underline{v}}_{n+\alpha_{\mathrm{m}}}}{\partial\underline{u}_{n+1}}\bm{\mathbf{M}}_{\alpha}+\frac{\partial\underline{u}_{n+\alpha_{\mathrm{f}}}}{\partial\underline{u}_{n+1}}\left(\bm{\mathbf{A}}^{\prime}(\underline{u}^{k}_{n+\alpha_{\mathrm{f}}})-\bm{\mathbf{M}}^{\prime}(\underline{u}^{k}_{n+\alpha_{\mathrm{f}}},\underline{p}_{\mathtt{c},{n+1}}^{k})\right)$
$\displaystyle=\rho_{0}\frac{\alpha_{\mathrm{m}}^{2}}{\alpha_{\mathrm{f}}\gamma^{2}\Delta
t^{2}}\bm{\mathbf{M}}_{\alpha}+\alpha_{\mathrm{f}}\left(\bm{\mathbf{A}}^{\prime}(\underline{u}^{k}_{n+\alpha_{\mathrm{f}}})-\bm{\mathbf{M}}^{\prime}(\underline{u}^{k}_{n+\alpha_{\mathrm{f}}},\underline{p}_{\mathtt{c},{n+1}}^{k})\right),$
(102)

with $\bm{\mathbf{A}}^{\prime}$, and $\bm{\mathbf{M}}^{\prime}$ being the
known tangent stiffness matrices from the quasi-stationary elasticity case,
see Eqs. 82 and 83. When using a coupling with the circulatory system we
compute the off diagonal matrices and lower right hand side, see Eqs. 84, 85
and 89, in terms of $\underline{u}^{k}_{n+\alpha_{\mathrm{f}}}$.

##  Appendix D Runge–Kutta–Fehlberg method

This methods is of order $\mathcal{O}(h^{4})$ but, by performing one extra
calculation, the error estimator is of order $\mathcal{O}(h^{5})$. This allows
for an adaptive stepsize that is determined automatically.

$\displaystyle y_{n+1}$ $\displaystyle=y_{n}+h\sum_{i=1}^{6}b_{1i}k_{i},$
$\displaystyle y^{*}_{n+1}$ $\displaystyle=y_{n}+h\sum_{i=1}^{6}b_{2i}k_{i},$
$\displaystyle e_{n+1}$
$\displaystyle=y_{n+1}-y^{*}_{n+1}=h\sum_{i=1}^{6}(b_{1i}-b_{2i})k_{i},$

where

$\displaystyle k_{1}$ $\displaystyle=f(t_{n},y_{n}),$ $\displaystyle k_{2}$
$\displaystyle=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1})),$ $\displaystyle k_{3}$
$\displaystyle=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2})),$
$\displaystyle\ \ \vdots$ $\displaystyle k_{6}$
$\displaystyle=f(t_{n}+c_{6}h,y_{n}+h(a_{61}k_{1}+a_{62}k_{2}+\cdots+a_{65}k_{5})).$

and coefficients given in Tab. 4.

$i$ | $c_{i}$ | $a_{i1}$ | $a_{i2}$ | $a_{i3}$ | $a_{i4}$ | $a_{i5}$ |
---|---|---|---|---|---|---|---
1 | 0 |  |  |  |  |  |
2 | $1/4$ | $1/4$ |  |  |  |  |
3 | $3/8$ | $3/32$ | $9/32$ |  |  |  |
4 | $12/13$ | $1932/2197$ | $-7200/2197$ | $7296/2197$ |  |  |
5 | $1$ | $439/216$ | $-8$ | $3680/513$ | $-845/4104$ |  |
6 | $1/2$ | $-8/27$ | $2$ | $-3544/2565$ | $1859/4104$ | $-11/40$ |
$j$ |  | $b_{j1}$ | $b_{j2}$ | $b_{j3}$ | $b_{j4}$ | $b_{j5}$ | $b_{j6}$
1 |  | $16/135$ | $0$ | $6656/12825$ | $28561/56430$ | $-9/50$ | $2/55$
2 |  | 25/216 | 0 | 1408/2565 | 2197/4104 | -1/5 | 0
Table 4: Butcher tableau for the Runge–Kutta–Fehlberg method.

##  Appendix E Direct Schur Complement Solver for a Small Number of
Constraints

Given the block system $\bm{\mathbf{A}}\in\mathbb{R}^{n\times n}$,
$\bm{\mathbf{D}}\in\mathbb{R}^{m\times m}$

$\displaystyle\begin{pmatrix}\bm{\mathbf{A}}&\bm{\mathbf{B}}\\\
\bm{\mathbf{C}}&\bm{\mathbf{D}}\end{pmatrix}\begin{pmatrix}\underline{x}\\\
\underline{y}\end{pmatrix}=-\begin{pmatrix}\underline{f}\\\
\underline{g}\end{pmatrix}$

with

$\displaystyle\bm{\mathbf{B}}=\begin{pmatrix}\underline{b}_{1}&\\!\\!\\!\vrule\\!\\!\\!&\cdots&\\!\\!\\!\vrule\\!\\!\\!&\underline{b}_{m}\end{pmatrix}\in\mathbb{R}^{n\times
m},\quad\bm{\mathbf{C}}=\begin{pmatrix}\underline{c}_{1}&\\!\\!\\!\vrule\\!\\!\\!&\cdots&\\!\\!\\!\vrule\\!\\!\\!&\underline{c}_{m}\end{pmatrix}^{\top}\in\mathbb{R}^{m\times
n},$

we can write the Schur complement system as

$\displaystyle(\bm{\mathbf{C}}\bm{\mathbf{A}}^{-1}\bm{\mathbf{B}}-\bm{\mathbf{D}})\underline{y}$
$\displaystyle=\underline{g}-\bm{\mathbf{C}}\bm{\mathbf{A}}^{-1}\underline{f}$
$\displaystyle\underline{x}$
$\displaystyle=\bm{\mathbf{A}}^{-1}\underline{f}-\bm{\mathbf{A}}^{-1}\bm{\mathbf{B}}\underline{y}.$

With

$\underline{r}=\bm{\mathbf{A}}^{-1}\underline{f},\quad\bm{\mathbf{S}}=\bm{\mathbf{A}}^{-1}\bm{\mathbf{B}}=\begin{pmatrix}\underline{s}_{1}&\\!\\!\\!\vrule\\!\\!\\!&\cdots&\\!\\!\\!\vrule\\!\\!\\!&\underline{s}_{m}\end{pmatrix}\in\mathbb{R}^{n\times
m},\ \underline{s}_{i}=\bm{\mathbf{A}}^{-1}\underline{b}_{i},\ i=1,\dots,m$
(103)

we get

$\displaystyle(\bm{\mathbf{C}}\bm{\mathbf{S}}-\bm{\mathbf{D}})\underline{y}$
$\displaystyle=\underline{g}-\bm{\mathbf{C}}\underline{r}$
$\displaystyle\underline{x}$
$\displaystyle=\underline{r}-\bm{\mathbf{S}}\underline{y}.$ (104)

The realization of (104) involves $m+1$ solves and the inversion of an
$m\times m$ matrix. Since $m$ is generally small this can be done
symbolically.

$\displaystyle{[\bm{\mathbf{C}}\bm{\mathbf{S}}]}_{ij}=\underline{c}_{i}\cdot\underline{s}_{j},\
\mbox{for}\ i,j=1,\dots,m.$

## References

  * Laslett et al. [2012] L. J. Laslett, P. Alagona, B. A. Clark, J. P. Drozda, F. Saldivar, S. R. Wilson, C. Poe, M. Hart,  The worldwide environment of cardiovascular disease: prevalence, diagnosis, therapy, and policy issues: a report from the american college of cardiology,  Journal of the American College of Cardiology 60 (2012) S1–S49.
  * Timmis et al. [2020] A. Timmis, N. Townsend, C. P. Gale, A. Torbica, M. Lettino, S. E. Petersen, E. A. Mossialos, A. P. Maggioni, D. Kazakiewicz, H. T. May, et al.,  European society of cardiology: cardiovascular disease statistics 2019,  European Heart Journal 41 (2020) 12–85.
  * Wilkins et al. [2017] E. Wilkins, L. Wilson, K. Wickramasinghe, P. Bhatnagar, J. Leal, R. Luengo-Fernandez, R. Burns, M. Rayner, N. Townsend,  European cardiovascular disease statistics 2017 (2017).
  * Smith et al. [2011] N. P. Smith, A. de Vecchi, M. McCormick, D. A. Nordsletten, O. Camara, a. F. Frangi, H. Delingette, M. Sermesant, J. Relan, N. Ayache, M. W. Krueger, W. H. W. Schulze, R. Hose, I. Valverde, P. Beerbaum, C. Staicu, M. Siebes, J. Spaan, P. J. Hunter, J. Weese, H. Lehmann, D. Chapelle, R. Rezavi,  euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling,  Interface Focus 1 (2011) 349–364.
  * Nordsletten et al. [2011] D. Nordsletten, M. Mccormick, P. J. Kilner, P. Hunter, D. Kay, N. P. Smith,  Fluid-solid coupling for the investigation of diastolic and systolic human left ventricular function,  International Journal for Numerical Methods in Biomedical Engineering 27 (2011) 1017–1039.
  * Karabelas et al. [2020] E. Karabelas, G. Haase, G. Plank, C. M. Augustin,  Versatile stabilized finite element formulations for nearly and fully incompressible solid mechanics,  Computational Mechanics 65 (2020) 193–215.
  * Elzinga and Westerhof [1973] G. Elzinga, N. Westerhof,  Pressure and flow generated by the left ventricle against different impedances,  Circulation Research 32 (1973) 178–186.
  * Liu et al. [2015] H. Liu, F. Liang, J. Wong, T. Fujiwara, W. Ye, K.-i. Tsubota, M. Sugawara,  Multi-scale modeling of hemodynamics in the cardiovascular system,  Acta Mechanica Sinica 31 (2015) 446–464.
  * Segers et al. [2008] P. Segers, E. Rietzschel, M. De Buyzere, N. Stergiopulos, N. Westerhof, L. Van Bortel, T. Gillebert, P. Verdonck,  Three-and four-element windkessel models: assessment of their fitting performance in a large cohort of healthy middle-aged individuals,  Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 222 (2008) 417–428.
  * Stergiopulos et al. [1999] N. Stergiopulos, B. E. Westerhof, N. Westerhof,  Total arterial inertance as the fourth element of the windkessel model,  American Journal of Physiology-Heart and Circulatory Physiology 276 (1999) H81–H88.
  * Wang et al. [2003] J.-J. Wang, A. B. O’Brien, N. G. Shrive, K. H. Parker, J. V. Tyberg,  Time-domain representation of ventricular-arterial coupling as a windkessel and wave system,  American Journal of Physiology-Heart and Circulatory Physiology 284 (2003) H1358–H1368.
  * Westerhof and Elzinga [1991] N. Westerhof, G. Elzinga,  Normalized input impedance and arterial decay time over heart period are independent of animal size,  American Journal of Physiology-Regulatory, Integrative and Comparative Physiology 261 (1991) R126–R133.
  * Westerhof and Stergiopulos [2000] N. Westerhof, N. Stergiopulos,  Models of the arterial tree.,  Studies in health technology and informatics 71 (2000) 65–77.
  * Alastruey et al. [2012] J. Alastruey, K. H. Parker, S. J. Sherwin,  Arterial pulse wave haemodynamics,  BHR Group - 11th International Conferences on Pressure Surges (2012) 401–442.
  * Blanco et al. [2014] P. J. Blanco, S. M. Watanabe, E. A. Dari, M. A. R. Passos, R. A. Feijóo,  Blood flow distribution in an anatomically detailed arterial network model: criteria and algorithms,  Biomechanics and Modeling in Mechanobiology 13 (2014) 1303–1330.
  * Formaggia et al. [2003] L. Formaggia, D. Lamponi, A. Quarteroni,  One-dimensional models for blood flow in arteries,  Journal of Engineering Mathematics 47 (2003) 251–276.
  * Mynard and Nithiarasu [2008] J. P. Mynard, P. Nithiarasu,  A 1D arterial blood flow model incorporating ventricular pressure, aortic vaive ana regional coronary flow using the locally conservative Galerkin (LCG) method,  Communications in Numerical Methods in Engineering (2008).
  * Müller and Toro [2014] L. O. Müller, E. F. Toro,  A global multiscale mathematical model for the human circulation with emphasis on the venous system,  International Journal for Numerical Methods in Biomedical Engineering 30 (2014) 681–725.
  * Marx et al. [2020] L. Marx, M. A. F. Gsell, A. Rund, F. Caforio, A. J. Prassl, G. Toth-Gayor, T. Kuehne, C. M. Augustin, G. Plank,  Personalization of electro-mechanical models of the pressure-overloaded left ventricle: fitting of Windkessel-type afterload models,  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 378 (2020) 20190342.
  * Bogaard et al. [2009] H. J. Bogaard, K. Abe, A. V. Noordegraaf, N. F. Voelkel,  The right ventricle under pressure: cellular and molecular mechanisms of right-heart failure in pulmonary hypertension,  Chest 135 (2009) 794–804.
  * Guazzi and Borlaug [2012] M. Guazzi, B. A. Borlaug,  Pulmonary hypertension due to left heart disease,  Circulation 126 (2012) 975–990.
  * Arts et al. [2005] T. Arts, T. Delhaas, P. Bovendeerd, X. Verbeek, F. Prinzen,  Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model,  Am. J. Physiol. Heart Circ. Physiol. 288 (2005) 1943–1954.
  * Blanco et al. [2010] P. J. Blanco, R. A. Feijóo, et al.,  A 3D-1D-0D Computational model for the entire cardiovascular system,  Computational Mechanics, eds. E. Dvorking, M. Goldschmit, M. Storti 29 (2010) 5887–5911.
  * Guidoboni et al. [2019] G. Guidoboni, L. Sala, M. Enayati, R. Sacco, M. Szopos, J. M. Keller, M. Popescu, L. Despins, V. H. Huxley, M. Skubic,  Cardiovascular function and ballistocardiogram: a relationship interpreted via mathematical modeling,  IEEE Transactions on Biomedical Engineering 66 (2019) 2906–2917.
  * Neal and Bassingthwaighte [2007] M. L. Neal, J. B. Bassingthwaighte,  Subject-specific model estimation of cardiac output and blood volume during hemorrhage,  Cardiovascular Engineering 7 (2007) 97–120.
  * Paeme et al. [2011] S. Paeme, K. T. Moorhead, J. G. Chase, B. Lambermont, P. Kolh, V. D’orio, L. Pierard, M. Moonen, P. Lancellotti, P. C. Dauby, et al.,  Mathematical multi-scale model of the cardiovascular system including mitral valve dynamics. application to ischemic mitral insufficiency,  Biomedical engineering online 10 (2011) 86.
  * Eriksson et al. [2013] T. S. E. Eriksson, A. J. Prassl, G. Plank, G. A. Holzapfel,  Influence of myocardial fiber/sheet orientations on left ventricular mechanical contraction,  Math Mech Solids 18 (2013) 592–606.
  * Kerckhoffs et al. [2007] R. C. Kerckhoffs, M. L. Neal, Q. Gu, J. B. Bassingthwaighte, J. H. Omens, A. D. McCulloch,  Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation,  Annals of biomedical engineering 35 (2007) 1–18.
  * Usyk et al. [2002] T. P. Usyk, I. J. LeGrice, A. D. McCulloch,  Computational model of three-dimensional cardiac electromechanics,  Computing and Visualization in Science 4 (2002) 249–257.
  * Fritz et al. [2014] T. Fritz, C. Wieners, G. Seemann, H. Steen, O. Dössel,  Simulation of the contraction of the ventricles in a human heart model including atria and pericardium,  Biomechanics and Modeling in Mechanobiology 13 (2014) 627–641.
  * Gurev et al. [2011] V. Gurev, T. Lee, J. Constantino, H. J. Arevalo, N. A. Trayanova,  Models of cardiac electromechanics based on individual hearts imaging data,  Biomechanics and Modeling in Mechanobiology 10 (2011) 295–306.
  * Gurev et al. [2015] V. Gurev, P. Pathmanathan, J.-L. Fattebert, H.-F. Wen, J. Magerlein, R. a. Gray, D. F. Richards, J. J. Rice,  A high-resolution computational model of the deforming human heart,  Biomechanics and Modeling in Mechanobiology 14 (2015) 829–849.
  * Hirschvogel et al. [2017] M. Hirschvogel, M. Bassilious, L. Jagschies, S. M. Wildhirt, M. W. Gee,  A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics,  International Journal for Numerical Methods in Biomedical Engineering 33 (2017).
  * Sainte-Marie et al. [2006] J. Sainte-Marie, D. Chapelle, R. Cimrman, M. Sorine,  Modeling and estimation of the cardiac electromechanical activity,  Computers and Structures 84 (2006) 1743–1759.
  * Augustin et al. [2016a] C. M. Augustin, A. Neic, M. Liebmann, A. J. Prassl, S. A. Niederer, G. Haase, G. Plank,  Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation,  Journal of Computational Physics 305 (2016a) 622–646.
  * Augustin et al. [2016b] C. M. Augustin, A. Crozier, A. Neic, A. J. Prassl, E. Karabelas, T. Ferreira da Silva, J. F. Fernandes, F. Campos, T. Kuehne, G. Plank, T. da Silva, J. F. Fernandes, F. Campos, T. Kuehne, G. Plank,  Patient-specific modeling of left ventricular electromechanics as a driver for haemodynamic analysis,  Europace 18 (2016b) iv121–iv129.
  * Crozier et al. [2016] A. Crozier, C. M. Augustin, A. Neic, A. J. Prassl, M. Holler, T. E. Fastl, A. Hennemuth, K. Bredies, T. Kuehne, M. J. Bishop, S. A. Niederer, G. Plank,  Image-Based Personalization of Cardiac Anatomy for Coupled Electromechanical Modeling,  Annals of Biomedical Engineering 44 (2016) 58–70.
  * Walmsley et al. [2015] J. Walmsley, T. Arts, N. Derval, P. Bordachar, H. Cochet, S. Ploux, F. W. Prinzen, T. Delhaas, J. Lumens,  Fast Simulation of Mechanical Heterogeneity in the Electrically Asynchronous Heart Using the MultiPatch Module,  PLOS Computational Biology 11 (2015) e1004284.
  * Prinzen et al. [2009] F. W. Prinzen, R. W. Mills, R. N. Cornelussen, L. J. Mulligan, M. Strik, L. M. Rademakers, N. D. Skadsberg, A. Van Hunnik, M. Kuiper, A. Lampert, T. Delhaas,  Left ventricular septal and left ventricular apical pacing chronically maintain Cardiac contractile coordination, pump function and efficiency,  Circulation: Arrhythmia and Electrophysiology 2 (2009) 571–579.
  * ver [2003] Intra-ventricular resynchronization for optimal left ventricular function during pacing in experimental left bundle branch block,  Journal of the American College of Cardiology 42 (2003) 558–567.
  * CIBC [2016] CIBC, 2016. URL: http://www.seg3d.org, seg3D: Volumetric Image Segmentation and Visualization. Scientific Computing and Imaging.
  * Neic et al. [2020] A. Neic, M. A. F. Gsell, E. Karabelas, A. J. Prassl, G. Plank,  Automating image-based mesh generation and manipulation tasks in cardiac modeling workflows using Meshtool,  SoftwareX 11 (2020) 100454.
  * Bayer et al. [2012] J. D. Bayer, R. C. Blake, G. Plank, N. A. Trayanova,  A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models,  Annals of Biomedical Engineering 40 (2012) 2243–2254.
  * Flory [1961] P. J. Flory,  Thermodynamic relations for high elastic materials,  Trans Faraday Soc 57 (1961) 829–838.
  * Land and Niederer [2018] S. Land, S. A. Niederer,  Influence of atrial contraction dynamics on cardiac function,  International journal for numerical methods in biomedical engineering 34 (2018) e2931.
  * Strocchi et al. [2020] M. Strocchi, M. A. Gsell, C. M. Augustin, O. Razeghi, C. H. Roney, A. J. Prassl, E. J. Vigmond, J. M. Behar, J. S. Gould, C. A. Rinaldi, M. J. Bishop, G. Plank, S. A. Niederer,  Simulating ventricular systolic motion in a four-chamber heart model with spatially varying robin boundary conditions to model the effect of the pericardium,  Journal of Biomechanics 101 (2020) 109645.
  * Usyk et al. [2000] T. P. Usyk, R. Mazhari, A. D. McCulloch,  Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle,  J. Elast. 61 (2000) 143–164.
  * Genet et al. [2014] M. Genet, L. C. Lee, R. Nguyen, H. Haraldsson, G. Acevedo-Bolton, Z. Zhang, L. Ge, K. Ordovas, S. Kozerke, J. M. Guccione,  Distribution of normal human left ventricular myofiber stress at end diastole and end systole: a target for in silico design of heart failure treatments,  Journal of Applied Physiology 117 (2014) 142–152.
  * Walker et al. [2005] J. C. Walker, M. B. Ratcliffe, P. Zhang, A. W. Wallace, B. Fata, E. W. Hsu, D. Saloner, J. M. Guccione,  MRI-based finite-element analysis of left ventricular aneurysm,  American Journal of Physiology-Heart and Circulatory Physiology 289 (2005) H692–H700.
  * Niederer et al. [2011] S. A. Niederer, G. Plank, P. Chinchapatnam, M. Ginks, P. Lamata, K. S. Rhode, C. A. Rinaldi, R. Razavi, N. P. Smith,  Length-dependent tension in the failing heart and the efficacy of cardiac resynchronization therapy,  Cardiovascular Research 89 (2011) 336.
  * Neic et al. [2017] A. Neic, F. O. Campos, A. J. Prassl, S. A. Niederer, M. J. Bishop, E. J. Vigmond, G. Plank,  Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model,  Journal of Computational Physics 346 (2017) 191–211.
  * Costa et al. [2013] C. M. Costa, E. Hoetzl, B. M. Rocha, A. J. Prassl, G. Plank,  Automatic Parameterization Strategy for Cardiac Electrophysiology Simulations.,  Computing in cardiology 40 (2013) 373–376.
  * ten Tusscher et al. [2004] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, A. V. Panfilov,  A model for human ventricular tissue.,  American Journal of Physiology. Heart and Circulatory Physiology 286 (2004) H1573–H1589.
  * Lumens et al. [2009] J. Lumens, T. Delhaas, B. Kirn, T. Arts,  Three-Wall Segment (TriSeg) Model Describing Mechanics and Hemodynamics of Ventricular Interaction,  Annals of Biomedical Engineering 37 (2009) 2234–2255.
  * Gurev et al. [2011] V. Gurev, T. Lee, J. Constantino, H. Arevalo, N. A. Trayanova,  Models of cardiac electromechanics based on individual hearts imaging data,  Biomechanics and modeling in mechanobiology 10 (2011) 295–306.
  * Rumpel and Schweizerhof [2003] T. Rumpel, K. Schweizerhof,  Volume-dependent pressure loading and its influence on the stability of structures,  International Journal for Numerical Methods in Engineering 56 (2003) 211–238.
  * Balay et al. [2018] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, a. Dave A. M, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, PETSc Users Manual, Technical Report ANL-95/11 - Revision 3.10, Argonne National Laboratory, 2018. URL: http://www.mcs.anl.gov/petsc.
  * Henson and Yang [2002] V. E. Henson, U. M. Yang,  BoomerAMG: A parallel algebraic multigrid solver and preconditioner,  in: Applied Numerical Mathematics, 2002\.
  * Vigmond et al. [2008] E. Vigmond, R. Weber dos Santos, A. Prassl, M. Deo, G. Plank,  Solvers for the cardiac bidomain equations,  Prog Biophys Mol Biol 96 (2008) 3–18.
  * Sugiura et al. [2012] S. Sugiura, T. Washio, A. Hatano, J. Okada, H. Watanabe, T. Hisada,  Multi-scale simulations of cardiac electrophysiology and mechanics using the University of Tokyo heart simulator,  Progress in Biophysics and Molecular Biology 110 (2012) 380–9.
  * Klotz et al. [2007] S. Klotz, M. L. Dickstein, D. Burkhoff,  A computational method of prediction of the end-diastolic pressure–volume relationship by single beat,  Nature Protocols 2 (2007) 2152–2158.
  * Arts et al. [1991] T. J. Arts, P. H. M. Bovendeerd, F. W. Prinzen, R. S. Reneman,  Relation between left ventricular cavity pressure and volume and systolic fibre stress and strain in the wall,  Biophys. J. 59 (1991) 93–102.
  * Firstenberg et al. [2001] M. S. Firstenberg, N. Greenberg, N. Smedira, P. McCarthy, M. J. Garcia, J. D. Thomas,  Noninvasive assessment of mitral inertness: clinical results with numerical model validation,  Comput. Cardiol. (2001) 613–616.
  * Otto [2006] C. M. Otto,  Valvular aortic stenosis: disease severity and timing of intervention.,  Journal of the American College of Cardiology 47 (2006) 2141–51.
  * Holzapfel [2000] G. A. Holzapfel, Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons Ltd, Chichester, 2000.
  * Deuflhard [2011] P. Deuflhard, Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms, volume 35, Springer Science & Business Media, 2011.
  * Press et al. [2007] W. H. Press, S. A. Teukolsky, W. T. Vettering, B. P. Flannery, Numerical Recipies: The Art of Scientific Computing, 3rd ed., Cambridge University Press, 2007\.
  * Kadapa et al. [2017] C. Kadapa, W. G. Dettmer, D. Perić,  On the advantages of using the first-order generalised-alpha scheme for structural dynamic problems,  Computers and Structures 193 (2017) 226–238.
  * Chung and Hulbert [1993] J. Chung, G. M. Hulbert,  A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalized-$\alpha$ Method,  Journal of Applied Mechanics 60 (1993) 371.
  * Newmark [1959] N. M. Newmark,  A method of computation for structural dynamics,  Journal of the engineering mechanics division 85 (1959) 67–94.