# Non-invasive estimation of left ventricle elastance using a multi-
compartment lumped parameter model and gradient-based optimization with
forward-mode automatic differentiation

Ryno Laubscher and Johan van der Merwe  Institute for Biomedical Engineering,
Department of Mechanical and Mechatronic Engineering, Stellenbosch University,
Stellenbosch, South Africa  Jacques Liebenberg and Philip Herbst  Division of
Cardiology, Faculty of Medicine and Health Sciences, Stellenbosch University,
Cape Town, South Africa

###### Abstract

Accurate estimates of left ventricle elastances based on non-invasive
measurements are required for clinical decision-making during treatment of
valvular diseases. The present study proposes a parameter discovery approach
based on a lumped parameter model of the cardiovascular system in conjunction
with optimization and non-invasive, clinical input measurements to approximate
important cardiac parameters, including left ventricle elastances. A subset of
parameters of a multi-compartment lumped parameter model was estimated using
1st order Adam gradient descent and hybrid Adam and 2nd order quasi-Newton
limited-memory Broyden-Fletcher-Goldfarb-Shanno optimization routines.
Forward-mode automatic differentiation was used to estimate the Jacobian
matrices and compared to the common finite differences approach. Synthetic
data of healthy and diseased hearts were generated as proxies for non-invasive
clinical measurements and used to evaluate the algorithm. Twelve parameters
including left ventricle elastances were selected for optimization based on
99% explained variation in mean left ventricle pressure and volume. The hybrid
optimization strategy yielded the best overall results compared to 1st order
optimization with automatic differentiation and finite difference approaches,
with mean absolute percentage errors ranging from 6.67% to 14,14%. Errors in
left ventricle elastance estimates for simulated aortic stenosis and mitral
regurgitation were smallest when including synthetic measurements for arterial
pressure and valvular flow rate at approximately 2% and degraded to roughly 5%
when including volume trends as well. However, the latter resulted in better
tracking of the left ventricle pressure waveforms and may be considered when
the necessary equipment is available.

###### keywords:

Cardiovascular simulation; Left ventricle elastances; Multi-compartment lumped
parameter model; Forward-mode automatic differentiation

††journal: Computers in Biology and Medicine

##  1 Introduction

The rapid rise of cardiovascular diseases (CVDs) in sub-Saharan Africa (SSA)
along with the high rates of communicable diseases in these regions, results
in a double burden experienced by the public health care sector, undermining
its effectiveness [1]. CVD related deaths in SSA account for approximately 11%
of all deaths in the region [2]. The most prevalent aetiologies of heart
failure in SSA that eventually result in morbidity and mortality are rheumatic
heart disease (RHD), ischemic heart disease and hypertensive heart disease
[3]. RHD is especially common in SSA as a result of low income status, poor
living conditions and limited access to medical care. The high prevalence and
late detection results in region-specific deaths that accounts for 23% of the
world’s deaths related to this disease. RHD usually involves the aortic or
mitral valve, resulting in either valvular stenosis or regurgitation. Long
standing neglected severe valve lesions ultimately leads to cardiac
decompensation and impairs the hearts ability to deliver oxygen rich blood to
various organs. Assessment of intrinsic ventricular function is of critical
importance in the management of patients with valve disease and is an
important trigger for recommending valve intervention, including valve
replacement surgery [4]. Unfortunately, the majority of clinically accessible
parameters used to assess ventricular function are measures of pump function
that are significantly load dependent, the most important of these, left
ventricular ejection fraction (LVEF). In the presence of severe mitral
regurgitation the intrinsic ventricular muscle function is difficult to judge
using LVEF as the LV offloads into a lower pressure left atrium (LA)
complicating its assessment. In severe aortic stenosis, intrinsic ventricular
muscle dysfunction from a number of associated comorbid conditions, can
similarly be difficult to differentiate from excessively high afterload
reducing muscle shortening. LV elastance, a load independent parameter, is an
appealing alternative and is considered the gold standard for LV systolic
function assessment and ventricular – arterial interaction [5]. To accurately
determine the LV elastance for a patient, invasive, multielectrode conductance
catheterisation is typically utilised [6]. In developing countries, this
invasive procedure is not readily available to clinicians, is expensive and is
not without risk. To non-invasively approximate a patient’s LV elastance,
clinicians typically use systolic and diastolic blood pressures, Doppler-based
measured stroke volume and derived LV ejection fraction with Chen’s formula
[7]. Using Chen’s formula with abnormal ventricular loading conditions, such
as those observed in patients with valvular diseases, leads to large
variability in the accuracy of the calculated LV elastance due to the load –
dependant nature of LVEF used in this calculation [8]. Therefore, a need still
exists to develop a robust approach to non-invasively approximate LV
elastances for diseased heart valves and ventricles.

An exciting, but experimental, approach to achieve robust non-invasive
estimation of patient LV elastance is the combination of cardiovascular lumped
parameter modelling (LPM), optimisation algorithms, and clinical measurements.
In this approach, the LPM acts as a mechanistic bridge between the non-
invasive measurements and the desired parameters (LV elastance). The
optimisation algorithms are used to minimise the difference between the LPM
results and the available measurements by tuning the mathematical model
parameters such as LV elastance, arterial resistance and pulmonary artery
compliance.

Recently, several researchers have investigated the feasibility of using LPMs
in conjunction with optimisation algorithms to find unknown cardiovascular
parameters. Herewith follows a brief overview of selected research papers on
this subject. Huang and Ying [9] developed an on-line parameter identification
algorithm to estimate the parameters of a 5 component arterial model. The
authors used synthetic data as a proxy for clinical measurements. The
synthetic data was generated using a multi-compartment LPM of a complete
cardiovascular system. To effectively approximate the parameters, the mean
squared error between the transient synthetic data and the corresponding model
predictions was minimized using a MATLAB Simulink model along with its
accompanying non-linear optimization package (fmincon [10]). The developed
algorithm was capable of recovering the unknown arterial resistances and
compliances from the synthetic data. Zhang et al. [11] developed a 0D-1D model
of a complete cardiovascular system and used clinically measured data to
identify key parameters such as arterial stiffnesses, peak flow time, distal
end scale factors and resistances. Non-invasive patient measurements were
performed to estimate artery diameters, ventricular stroke volumes, blood
pressure waveforms and axial blood flow velocities. To find the unknown
parameter set, the mean squared differences between the model predictions and
measured brachial arterial pressures and aortic velocities were minimized
using the Levenberg-Marquardt optimization algorithm. Seeing as only non-
invasive stroke volumes and brachial pressures were used by the authors, it
limits the approach to persons with no heart diseases such as aortic valve
stenosis. It should be stated that modelling heart diseases was not the
purpose of the previous authors’ research. Colunga et al. [12] used invasive
right heart catheterisation data and LPMs to identify cardiovascular
parameters such as pulmonary and systemic resistances, pulmonary valve
resistance, pulmonary and systemic arterial compliances and right ventricle
elastances. The LPM developed by the authors consisted of 6 components (e.g.
systemic artery, left ventricle, pulmonary artery) and the squared differences
between the predicted and clinically measured arterial blood pressures, right
heart ventricle pressures and cardiac outputs were minimised using the
Levenberg-Marquardt optimisation routine in MATLAB for various patients. The
optimised LPMs could accurately recreate measured waveforms but seeing as the
identified parameters were not directly measured, the absolute accuracy of the
identified parameters is unknown. Keshavarz-Motamed [13] developed a
comprehensive workflow coined C3VI-CMF, to identify circulatory system
parameters and heart functions using MATLAB Simscape and the non-linear design
optimisation toolbox in MATLAB Simulink. The developed tool uses Doppler
echography and arm cuff pressure measurements along with a multi-compartment
LPM to identify aortic compliance, systemic arterial and venous compliance,
systemic arterial resistance and mean pulmonary valve flow rate. The remainder
of the LPM parameters, such as ventricular elastances, were set to typical
clinical values. The ventricular pressure waveforms predicted by the LPM with
the identified parameters were compared to the ventricular catheterisation
data and showed excellent agreement. In [13], errors between the predicted LPM
outputs and the non-invasive measurements were optimised using the fmincon
function in MATLAB. Bjørdalsbakke et al. [14] developed a program capable of
estimating systemic loop parameters such as systemic compliances, LV
elastances (diastolic and systolic), systemic resistance and aorta blood
inertia using the trust region reflective algorithm [15] from the free and
opensource Python library SciPy. The developed LPM did not consider the
pulmonary loop and approximated the systemic loop using a 3-element component.
The authors developed an innovative technique called the stepwise subset
reduction method to reduce the effect of less important parameters during
optimisation. Synthetic data was generated using the LPM model and the
difference between the model predictions and the synthetic data was used to
identify the unknown parameters. The achieved mean absolute percentage errors
(MAPEs) ranged between 1-10% for the identified parameters.

For the previously discussed MATLAB- and Python-based models, the parameter-
cost function Jacobian matrices used in the optimisation algorithms are
estimated using the finite differences approach [10], [16]. Gradients obtained
using finite differences are typically computationally expensive for large
problems, such as system-level transient simulations, and tend to be
numerically unstable seeing as it only approximates the gradients [17]. An
alternative approach to estimate the Jacobian matrices is to use a method
typically applied in deep learning called automatic differentiation (autodiff)
[18], which calculates the analytical gradients using the chain rule and
computational graphs constructed from the mathematical operations in the
computer model. In the present work, we propose a parameter discovery
algorithm which leverages automatic differentiation through the unsteady
numerical simulation procedure of a cardiovascular LPM to identify important
parameters, relating to and including LV elastances for diseased and healthy
left heart valves. The algorithm consists of a multi-compartment LPM of an
entire cardiovascular system and 1st and 2nd order gradient-descent
optimisation routines. To demonstrate the performance and accuracy of the
proposed algorithm, synthetic data generated using the LPM is used as a proxy
for clinical measurements. Only measurements which can be non-invasively
acquired such as ventricular volumes, valvular flow rates and arterial
pressures are utilised to discover the unknown system parameters. To identify
a parameter subset to be optimised, a local sensitivity study is performed
also using automatic differentiation to calculate the LV pressure and volume
sensitivities relative to the cardiovascular model parameters. To lower costs
of deploying the proposed algorithm and to enable reproducibility, the
computer models are developed in free and opensource Julia 1.7.0 libraries
namely DifferentialEquations.jl [19], ForwardDiff.jl [20] (automatic
differentiation), Optim.jl [21] (optimisation framework) and Flux.jl [22]
(first-order gradient-descent optimisers).

##  2 Materials and methods

###  2.1 Overview of proposed model

In the present work, synthetic data of ventricular volume changes, valvular
flow rates and arterial pressure waveforms for a single beat are generated
using a multi-compartment LPM and serves as substitutions for actual clinical
measurements. The benefit of using synthetic data is that the true parameters
used to generate the measurements are known, thus the accuracy of the
parameter optimisation model can be evaluated. The objective of the proposed
model is to recover the LV elastances used during synthetic data generation,
by utilising the substitute non-invasive measurements and the LPM equations.
The development of the proposed model can be divided into two parts namely
parameter subset selection and parameter optimisation (figure 1).

Figure 1: Computer model flow chart

Seeing as various cardiovascular parameters such as arterial resistances and
unstressed LV volume influence the pressure-volume relationship of the LV and
is typically unknown prior to invasive measurements, a sensitivity analysis
using the LPM is performed. This is done to define a reduced parameter subset
$\bar{\theta}$ which contains the LV elastances and important parameters which
influence the LV pressure-volume loop. In the current work, the remainder of
the model parameters which is not to be optimised are fixed to typical values
found in literature. More on this later on in the paper. Once the parameter
subset has been identified, it can be optimised by minimizing the error
between the parameter-specific model predictions and the synthetic data.

In figure 1, the program flow chart for the parameter optimisation is
depicted. Initially $\bar{\theta}$ is randomly initialised to values within
typical clinical ranges. Next, the parameters are used along with the LPM and
an ODE solver to produce a solution for a single heartbeat period
$\bar{X}=M(\bar{\theta})$. Using the parameter-specific mathematical model
solution, the relevant waveforms $\tilde{X}$ is extracted from $\bar{X}$ and
along with the synthetic data $\hat{X}$ fed into a cost function $C$. If the
resulting cost function value is below the user specified convergence criteria
$\epsilon$ the program exits. Conversely, if the cost function value is above
$\epsilon$, the gradients of the cost function with respect to the model
parameters $\frac{\partial C(\hat{X},\bar{X})}{\partial\bar{\theta}}$ is
calculated using forward-mode automatic differentiation. Next, $\frac{\partial
C}{\partial\bar{\theta}}$ is used along with a gradient-descent optimiser to
calculate the updated parameters $\bar{\theta}^{\text{new}}$ new which is then
fed back to the LPM and ODE solver and the process is repeated.

In the following subsections, the lumped parameter modelling, parameter
optimisation and sensitivity analysis (reduced parameter subset selection)
methods will be discussed in more detail.

###  2.2 Mathematical modelling and synthetic datasets

####  2.2.1 Lumped parameter modelling

The basis of the present work is a multi-compartment LPM of the entire
cardiovascular system, which is depicted in figure 2 below. The LPM developed
in the present work consists of the four heart chambers and their respective
downstream valves, the systemic loop and the pulmonary loop.

Figure 2: Cardiovascular multi-compartment lumped parameter model. Aortic –
AO, mitral – MI, pulmonary – PO, tricuspid – TI, aortic sinus – AS, systemic
arteries – SAT, systemic arterioles – SAR, systemic capillaries – SCP,
systemic veins – SVN, pulmonary sinus – PS, pulmonary arteries – PAT,
pulmonary arterioles - PAR, pulmonary capillaries – PCP, pulmonary veins - PVN

The pressure-volume relationship of the heart chambers can be mathematically
modelled using the Suga et al. [23] relation as seen in equation 1. In
equation 1, $P_{lv}(t)$ [mmHg] is the LV pressure at time $t$, $P_{lv,0}$ is
the unstressed LV pressure (set to a value of 1 for all heart chambers [24]),
$e_{lv}(t)$ [s] is the LV time-varying elastance function, $V_{lv}(t)$ [mL] is
the instantaneous LV volume and $V_{lv,0}$ is the unstressed LV volume.

$P_{lv}(t)=P_{lv,0}+e_{lv}(t)\left(V_{lv}(t)-V_{lv,0}\right)$ (1)

The ventricle volume can be simulated using a simple mass balance relation
over the ventricle control volume as shown in equation 2. In equation 2,
$Q_{mi}(t)$ [mL/s] is the volume flow rate through the mitral valve and
$Q_{ao}(t)$ is the flow rate through the aortic valve.

$\frac{dV_{lv}}{dt}=Q_{mi}(t)-Q_{ao}(t)$ (2)

The time-varying elastance function for the LV is calculated using the
systolic and diastolic elastances $E_{lv,s}$, $E_{lv,d}$ [mmHg/mL] and the
contractility activation function $f(t)$ as shown in equation 3.

$e_{lv}(t)=E_{lv,d}+\frac{E_{lv,s}-E_{lv,d}}{2}\cdot f(t)$ (3)

In the present work the ventricle activation function used was taken from the
work of Bozkurt [25] and is shown in equation 4.

$f_{i}(t)=\begin{cases}1-\cos\left[\left(\frac{t}{T_{1}}\pi\right)\right],&\text{if
}0\leq t<T_{1}\\\ 1+\cos\left[\frac{t-T_{1}}{T_{2}-T_{1}}\pi\right],&\text{if
}T_{1}\leq t<T_{2}\\\ 0,&\text{if }T_{2}\leq t<T\par\end{cases}$ (4)

In equation 4, $T_{1}=0.3T$ [s] is the time at end of systole, $T_{2}=0.45T$
is the time at end of ventricular relaxation [25] and $T$ is the heart beat
period which in the present work was fixed to 1 [s]. These time intervals can
typically be acquired from electrocardiogram (ECG) and Doppler
echocardiography (DE) data. To model the right ventricle, equations 1-6 are
similarly applied but with corresponding right heart parameters such as
$E_{rv,s}$ and $E_{rv,d}$ and tricuspid and pulmonary valve flow rates
$Q_{ti}(t)$ and $Q_{po}(t)$.

The atria pressures and instantaneous chamber volumes are calculated similarly
to the ventricles (equations 1 and 2), but with different definitions for
time-varying elastance functions and contractility activation functions. The
time-varying left atrial (LA) elastance function is calculated using equation
5. The right atrial elastance function is similarly calculated.

$e_{la}(t)=E_{la,min}+\frac{E_{la,max}-E_{la,min}}{2}\cdot f_{a}(t-D)$ (5)

In equation 5, $E_{la,min}$ and $E_{la,max}$ are the minimal and maximal LA
elastances, $f_{a}(t)$ is the atrial contractility activation function which
in turn is calculated using equation 6 and $D=0.04$ [s] is the time of atrial
relaxation.

$f_{a}(t)=\begin{cases}0,&\text{if }0\leq t<T_{a}\\\
1-\cos\left[2\pi\frac{t-T_{a}}{T-T_{a}}\right],&\text{if }T_{a}\leq
t<T\end{cases}$ (6)

In equation 6, $T_{a}=0.8T$ is the time at onset of atrial contraction.
Similar to the ventricular time intervals, $T_{a}$ and $D$ can be measured
using ECG and DE data. Table 1 below contains the nominal parameters [26] used
in the LPM.

Table 1: Heart model nominal parameters (values in parenthesis indicate upper and lower boundaries for sensitivity analysis and normalisation) Parameters | Left heart | Right heart
---|---|---
| _Atria_ |
$E_{max}$ $\left[\frac{mmHg}{mL}\right]$ | $0.25(0.0,1.0)$ | $0.15(0.0,1.0)$
$E_{min}$ $\left[\frac{mmHg}{mL}\right]$ | $0.15(0.0,0.5)$ | $0.15(0.0,0.5)$
$V_{0}$ $\left[mL\right]$ | $4(1.0,20.0)$ | $4(1.0,20.0)$
| _Ventricles_ |
$E_{d}$ $\left[\frac{mmHg}{mL}\right]$ | $0.1(0.0,1.0)$ | $0.1(0.0,0.5)$
$E_{s}$ $\left[\frac{mmHg}{mL}\right]$ | $2.5(0.5,5.0)$ | $1.15(0.5,5.0)$
$V_{0}$ $\left[mL\right]$ | $5(1.0,20.0)$ | $10(1.0,50.0)$

The pressure-volume flow rate relation of the heart valve models was simulated
using a orifice model proposed by Korakianitis and Shi [24]. For this heart
valve model the flow rate through the $i^{t}h$ valve is calculated using
equation 7, where $i=ao,mi,ti,po$ (aortic, mitral, tricuspid and pulmonary).

$Q_{i}(t)=CQ\cdot A_{r}(t)\sqrt{\Delta P(t)}$ (7)

In equation 7, $CQ$ is the flow coefficient which for the atrioventricular
valves is set to 400 and for the semilunar valves to 350, $A_{r}(t)$ is the
time-dependent area opening fraction of the valve, and $\Delta P(t)$ the
pressure drop across the valve which in turn is calculated using equation 8.

$\Delta P(t)=\begin{cases}P_{in}(t)-P_{ex}(t),\text{if }&P_{in}\geq P_{ex}\\\
P_{ex}(t)-P_{in}(t),\text{if }&P_{in}<P_{ex}\end{cases}$ (8)

In equation 8, $P_{in}(t)$ is the instantaneous valve inlet pressure and
$P_{ex}(t)$ the exit or downstream pressure. For example, the aortic valve
inlet pressure would be the LV pressure $P_{lv}(t$ and the exit pressure the
aortic sinus pressure $P_{AS}(t)$ as shown in figure 2. The valve area opening
fraction is simulated using simple diode behaviour, where $A_{r}$ either takes
a value of 0 or 1 depending on the pressure gradient across the specific heart
valve as shown in equation 9.

$A_{r}(t)=\begin{cases}1,\text{if }&P_{in}(t)\geq P_{ex}(t)\\\ 0,\text{if
}&P_{in}(t)<P_{ex}(t)\end{cases}$ (9)

The systemic and pulmonary loops each consists of five components which are
sinuses, arteries, arterioles, capillaries and veins. For each of these
compartments a set of equations are used to simulate their internal pressure
and flow rate dynamics. Note that only the systemic loop equations will be
discussed below, but the pulmonary loops solves similar equations with
different parameters. The flow rate through the aortic sinus and its
corresponding inlet pressure are simulated using two ordinary differential
equations as shown as follows in equations 10 and 11.

$L_{AS}\frac{dQ_{AS}}{dt}=\left(P_{AS}-P_{SAT}\right)-R_{AS}Q_{AS}$ (10)

$C_{AS}\frac{dP_{AS}}{dt}=Q_{AO}-Q_{AS}$ (11)

In equations 10 and 11, $L_{AS}$ $[\frac{mmHg\cdot s^{2}}{mL}]$ is the blood
flow inertia through the aortic sinus, $Q_{AS}$ is the flow rate of blood
through the sinus, $P_{AS}$ is the inlet sinus pressure, $P_{SAT}$ is the
systemic arterial inlet pressure, $R_{AS}$ $[\frac{mmHg\cdot s}{mL}]$ is the
sinus flow resistance and $C_{AS}$ $[\frac{mL}{mmHg}]$ is the sinus vessel
wall compliance. The systemic arteries are simulated similarly to the aortic
sinus. The equations governing the arterial pressure and flow rate are given
in equations 12 and 13.

$L_{SAT}\frac{dQ_{SAT}}{dt}=\left(P_{SAT}-P_{SAR}\right)-R_{SAT}Q_{SAT}$ (12)

$C_{SAT}\frac{dP_{SAT}}{dt}=Q_{AS}-Q_{SAT}$ (13)

The systemic arterioles and capillaries are simulated using only hydraulic
resistance components due to the rigidity and small diameter of the vessel
walls [27]. The flow rates through these sections of the cardiovascular system
are assumed to be steady, thus due to the mass conservation law the flow rates
are all assumed to be equal $Q_{SAT}=Q_{SAR}=Q_{SCP}$. The pressure drop
through these sections can, therefore, be calculated as:

$Q_{SAT}R_{SAR}=P_{SAR}-P_{SCP}$ (14)

$Q_{SAT}R_{SCP}=P_{SCP}-P_{SCN}$ (15)

The systemic veins are simulated using compliance-resistance components with
negligible blood inertia [24]. The governing differential equation for the
systemic vein inlet pressure can be seen in equation 16.

$C_{SVN}\frac{dP_{SVN}}{dt}=Q_{SCP}-Q_{SVN}$ (16)

To calculate the flow rate through the venous system the hydraulic resistance
expression is utilised. The driving pressure for the systemic vein flow rates
is the difference in systemic vein inlet pressure and the right atrium
pressure $P_{ra}(t)$, as seen in equation 17.

$Q_{SVN}R_{SVN}=P_{SVN}-P_{ra}(t)$ (17)

The various model parameters used in the previous equations (10-17) for both
the systemic and pulmonary loops can be seen in table 2.

Table 2: Systemic and pulmonary loop parameters (values in parenthesis indicate upper and lower boundaries for sensitivity analysis and normalisation) [24] Components |  $R$ $\left[\frac{mmHg\cdot s}{mL}\right]$ |  $L$ $\left[\frac{mmHg\cdot s^{2}}{mL}\right]$ |  $C$ $\left[\frac{mL}{mmHg}\right]$
---|---|---|---
_Systemic network_ |  |  |
AS | 0.003 ($3\cdot 10^{-4}$, 0.03) |  $6.2\cdot 10^{-5}$ ($1\cdot 10^{-4}$, 0.001) | 0.08 (0.008,0.8)
SAT | 0.05 (0.005,1.0) | 0.0017 ($1.7\cdot 10^{-3}$, 0.017) | 1.6 (0.16,3.2)
SAR | 0.5 (0.05, 1.0) | - | -
SCP | 0.52 (0.05, 1.0) | - | -
SVN | 0.075 (0.0075, 0.75) | - | 20.5 (5.0, 50.0)
_Pulmonary network_ |  |  |
PS | 0.002 (0.0002, 0.002) |  $5.2\cdot 10^{-5}$ ($1\cdot 10^{-4}$, $1\cdot 10^{-3}$) | 0.18 (0.018, 2.0)
PAT | 0.01 (0.001, 0.1) | 0.0017 ($1.7\cdot 10^{-3}$, 0.017) | 3.8 (0.038, 6.0)
PAR | 0.05 (0.005, 1.0) | - | -
PCP | 0.05 (0.025, 2.0) | - | -
PVN | 0.006 ($6\cdot 10^{-4}$, 0.01) | - | 20.5 (5.0, 50.0)

To simulate the dynamics of the cardiovascular system, the differential-
algebraic equations (DAEs) for the time-varying elastance model (equations
1-6), the heart valve model equations (equations 7-9) and vasculature
equations (equations 10-17) for the left and right heart chambers, valves and
systemic and pulmonary loops are solved simultaneously using an ODE
integrator. For the current work the Tsitouras 5/4 Runge-Kutta variable time-
step method [28] is applied with relative and absolute tolerances set to 1E-4
and 1E-6 respectively. In total the model solves for 14 dependent variables
and uses the parameters from tables 1 and 2 as model inputs along with the
initial conditions vector $\bar{X}_{init}$. These 14 dependent variables are
defined by the model equations for the four heart chambers, four heart valves
and systemic and pulmonary loops as discussed in section 2.2.1. For a more
detailed discussion of the model equations the interested reader is referred
to [29]. Once the model is solved, the corresponding output has the size
$M(\tilde{\theta},\bar{X}_{i}nit)=\tilde{X}\in\mathcal{R}^{14\times N_{t}}$
where $\tilde{\theta}$ is a vector containing all the parameters mentioned
above, $\tilde{X}$ is the model solution matrix and $N_{t}$ is the total
number of time steps taken by the ODE integrator over the simulation time
$T=1$ [s]. Note, that the parameters to be optimised is only a subset of the
complete parameter vector thus $\bar{\theta}\in\bar{\theta}$ and similarly
$\bar{X}$ is only a subset of the solution matrix thus $\bar{X}\in\tilde{X}$.

####  2.2.2 Synthetic data and lumped parameter model initial conditions

Clinical measurement data is typically discrete values sampled at a certain
frequency, therefore, the synthetic data matrix has a size of
$\hat{X}\in\mathcal{R}^{d\times N_{t}^{*}}$, where $d$ is the number of
different types of measurements and $N_{t}^{*}$ the number of sampled time
steps for each performed measurement. In the present work $N_{t}^{*}=100$
which equivalent to a sampling frequency of 100 Hz. To compare the simulated
trends to the synthetic data, a subset of the simulated values is sampled from
$\tilde{X}$ at time steps corresponding to the synthetically measured time
stamps, thus resulting in the simulation subset
$\bar{X}\in\mathcal{R}^{d\times N_{t}^{*}}$.

In the current work, three measurement datasets containing different
combinations of typical clinical non-invasive measurements are evaluated to
find the best-performing configuration that results in the smallest error
between the actual model parameters and optimised ones. Table 3 shows the
dimensions and measurement collections of these three $\hat{X}$ datasets.

Table 3: Time-varying elastance model parameters for atriums. Ref. [25], [26] Description | Dimensions | Simulated measurements
---|---|---
Dataset 1 (D1) |  $d=5$, $N_{t}^{*}=100$ |  $P_{SAT}$, $V_{la}$, $V_{lv}$, $V_{ra}$, $V_{rv}$
Dataset 2 (D2) |  $d=5$, $N_{t}^{*}=100$ |  $P_{SAT}$, $Q_{ao}$, $Q_{mi}$, $Q_{po}$, $Q_{ti}$
Dataset 3 (D3) |  $d=9$, $N_{t}^{*}=100$ | D1 + D2

To measure the clinical trends in table 3, different medical equipment can be
used. For the current study the following measurements are proposed for
further clinical studies. For the heart chamber volume trends, 3D echography
or magnetic resonance imaging (MRI) could be utilised and Doppler echography
to measure the heart valve flow rates. The brachial arterial pressure can be
continuously measured using a CNAP Monitor [11] and the volume-clamp method.

To simulate the dynamics of the LPM dependent variables such as LV volume
changes, systemic arterial pressures and venous pressures, using the model
equations, 14 initial conditions are required. The initial conditions vector
can be defined as
$\bar{X}_{init}=[V_{lv}^{init},V_{la}^{init},P_{AS}^{init},Q_{AS}^{init},P_{SAT}^{init},$

$Q_{SAT}^{init},P_{SVN}^{init},V_{rv}^{init},V_{ra}^{init},P_{PS}^{init},Q_{PS}^{init},P_{PAT}^{init},Q_{PAT}^{init},P_{PVN}^{init}]$.
In a clinical application of the current work, these initial conditions should
also be extracted from non-invasive measurements. The initial heart chamber
volumes, $V_{lv}^{init},V_{la}^{init},V_{rv}^{init}$ and $V_{ra}^{init}$, can
be directly taken as the first entries in the corresponding synthetic data
trends. The initial aortic and pulmonary sinus flow rates, $Q_{AS}^{init}$ and
$Q_{PS}^{init}$ , are assumed to be equal to the initial entries in the
synthetic datasets for the aortic and pulmonary valvular flow rates. The
initial systemic arterial pressure $P_{SAT}^{init}$ can be extracted from the
synthetic dataset trends. The initial aortic sinus pressure $P_{AS}^{init}$ is
assumed to be equal to the initial systemic arterial pressure. For the
remaining initial conditions, additional clinical measurements are required.
The total systemic and pulmonary arterial flow rates $Q_{SAT}^{init}$ and
$Q_{PAT}^{init}$ is estimated using the LV and RV stroke volume measured using
Doppler echography and divided by the heartbeat interval time. The initial
pulmonary and systemic venous pressures $P_{SVN}^{init}$ and $P_{PVN}^{init}$
are important initial conditions, seeing as it drives the blood flow into the
atria. These pressures are difficult to estimate non-invasively, therefore, in
the present work the initial venous pressures are optimised along with the
other model parameters. Similarly, the pulmonary sinus and arterial pressure
initial conditions $P_{PS}^{init}$ and $P_{PAT}^{init}$ are also difficult to
estimate non-invasively, and therefore, are also optimised along with the
model parameters.

###  2.3 Parameter optimisation

The parameter vector containing the unknown cardiovascular system quantities
is in essence a row vector of values which individually each have a
significant influence on the LV pressure-volume relationship. How the
parameters in this vector are selected will be discussed in section 2.4. An
example of the parameter vector is shown in equation 18.

$\bar{\theta}=[E_{lv,s},E_{lv,d},…,R_{AS},…,C_{PVN}]$ (18)

To find $\bar{\theta}$ which minimizes the differences between the simulation
results and measurement data, a sum squared error (SSE) cost function is
defined for each of the measurements, as shown in equation 19 for the
$j=1,2,..,d$ measurement types (e.g. LV volume, arterial pressure or aortic
valve flow rate).

$J(\bar{x}_{j},\hat{x}_{j})=\left(\sum_{i=1}^{N_{t}^{*}}(\bar{x}_{j}^{i}(\bar{\theta})-\hat{x}_{j}^{i})^{2}\right)_{j}$
(19)

In equation 19, $\bar{x_{j}^{i}}$ is the $j^{th}$ simulation output at time
step $i$ and $\hat{x}_{j}^{i}$ is the $j^{th}$ synthetically measured value
(e.g. arterial pressure or LV volume) at time step $i$. The objective of the
optimisation routine is to minimise the sum of the individual measurement type
cost functions as shown in equation 20.

$C(\bar{X}(\bar{\theta}),\hat{X})=\sum_{i=1}^{d}\left(J(\bar{x}_{i}(\bar{\theta}),\hat{x}_{i})\right)$
(20)

To speed up model convergence, the parameter vector entries are normalised
using min-max scaling and user specified upper and lower parameter boundaries
(see tables 1 and 2). Min-max scaling results in parameters having numeric
values between 0-1. The scaling transformation can be seen in equation 21,
where $\bar{\theta}^{*}$ is the scaled parameter vector, $\bar{\theta}_{lb}$
is a vector of lower boundary parameter values and $\bar{\theta}_{ub}$ is a
vector of upper boundary parameter values. Note that during the optimisation
routine, the parameter vector is rescaled before being sent to the LPM and ODE
integrator routine.

$\bar{\theta}^{*}=\frac{\bar{\theta}-\bar{\theta}_{lb}}{\bar{\theta}_{ub}-\bar{\theta}_{lb}}$
(21)

Due to the difference in magnitudes of volumes, flow rates and pressures, the
simulation outputs and measurement data are also scaled to values between 0-1
before being used in the cost function calculations. Min-max scaling is again
used along with the maximum and minimum measured values for each measurement
type. For example, the scaled simulated arterial pressures can be calculated
using equation 22, where the subscript $p$ indicates arterial pressure values
being compared from the simulation output and synthetic measurement vectors
$\bar{x}$ and $\hat{x}$ respectively.

$\bar{x}_{p}^{*}=\frac{\bar{x}_{p}-\text{min}(\bar{x_{p}})}{\text{max}(\bar{x_{p}})-\text{min}(\bar{x_{p}})}$
(22)

To find the parameters which minimises the cost function $C(\bar{\theta})$, a
combination of first- and second-order gradient descent optimisers are
utilised. The adaptive moment estimation (Adam) first-order optimiser [30] was
selected. The Adam optimiser algorithm is shown in equation 23.

$\begin{gathered}\bar{m}\leftarrow\beta_{1}\bar{m}+(1-\beta_{1})\nabla_{\theta}C(\bar{\theta})\\\
\bar{s}\leftarrow\beta_{1}\bar{s}+(1-\beta_{2})\nabla_{\theta}C(\bar{\theta})\otimes\nabla_{\theta}C(\bar{\theta})\\\
\bar{m}\leftarrow\frac{\bar{m}}{1-\beta_{1}^{t}}\\\
\bar{s}\leftarrow\frac{\bar{s}}{1-\beta_{2}^{t}}\\\
\bar{\theta}^{\text{new}}\leftarrow\bar{\theta}-\eta\bar{m}\otimes\sqrt{(\bar{s}+\epsilon)^{-1}}\end{gathered}$
(23)

The scaling $\bar{s}$ and momentum $\bar{m}$ matrices are initialized to 0 at
the start of the Adam training algorithm, $t$ is the iteration counter,
$\epsilon=10^{-8}$ is the smoothing term, $\beta_{1}$ is the momentum decay
hyperparameter and is set to 0.9 and $\beta_{2}$ is the scaling hyperparameter
and is set to 0.999. In equation 23, $\nabla_{\theta}C(\bar{\theta})$ are the
gradients of the cost function with respect to the optimisation parameters.
For the optimisation runs, the learning rate parameter $\eta$ is fixed to a
value of 0.005. In the development of the optimisation model, other optimisers
were tested such as the standard gradient descent, momentum and Nestorov
algorithms, but to stabilise these simulations, a small learning rate of 1E-5
had to be used. Therefore, it was decided to only use the Adam optimiser due
to its superior stability for the current problem.

As mentioned, in addition to using first-order optimisation, a second-order
optimiser was also applied to find the optimal parameters. Second-order
optimisers use the Hessian in addition to the gradients to determine the
parameter step direction and step lengths during optimisation [18]. In the
current work, the quasi-Newton Limited-memory Broyden-Fletcher-Goldfarb-Shanno
(L-BFGS) [31] optimiser was utilised, which does not directly calculate the
Hessian but rather approximates it using the first-order gradients.

For both the above mentioned first- and second-order optimisers, the cost
function gradients with respect to the parameters
$\nabla_{\theta}C(\bar{\theta})$, are required. Forward-mode automatic
differentiation is applied to calculate these gradients directly by tracking
both values and derivatives as they are passed through the computational
graphs created by the various mathematical operations. The ForwardDiff.jl
library used in the present work, is capable of traversing any native Julia
code, therefore, it is capable of differentiating through the ODE integrator
to calculate the required gradients in a computationally efficient manner, as
will be demonstrated in the results section.

###  2.4 Sensitivity analysis

To identify the model parameters (tables 1 and 2) which have a significant
influence on the LV behaviour, a local sensitivity analysis is performed,
excluding model initial conditions. This entails calculating the relative
sensitivity percentages for each parameter using the LPM, the nominal values,
upper and lower parameter boundary values and forward-mode automatic
differentiation through the LPM solution.

The main outputs which underpin the LV contractility is its time-dependent
volume and internal pressure changes. Therefore, to calculate the parameter
sensitivity of these values, the gradients of the mean pressure and volume
with respect to each parameter, listed in tables 1 and 2, are calculated
centred around the nominal values. Seeing that the model parameters vary
significantly in order of magnitude and units, the mean LV pressure and volume
parameter gradients were multiplied by their respective parameter ranges as
seen below. Equations 24 and 25, shows the mean LV pressure and volume
sensitivity indices calculated for the $i^{th}$ parameter $\theta_{i}$.

$SI_{P,i}=\left|\frac{\partial\left(\frac{1}{N_{t}^{*}}\sum_{j=1}^{N_{t}^{*}}P_{lv}^{j}\right)}{\partial\theta_{i}}\right|\cdot\left(\theta_{i,ub}-\theta_{i,lb}\right)$
(24)

$SI_{V,i}=\left|\frac{\partial\left(\frac{1}{N_{t}^{*}}\sum_{j=1}^{N_{t}^{*}}V_{lv}^{j}\right)}{\partial\theta_{i}}\right|\cdot\left(\theta_{i,ub}-\theta_{i,lb}\right)$
(25)

In equations 24 and 25, $P_{lv}^{j}$ is the LPM simulated LV pressure at time
sample $j$, similarly $V_{lv}^{j}$ is the simulated LV volume, $\theta_{i,ub}$
is the upper parameter boundary value for parameter $i$ and similarly
$\theta_{i,lb}$ the lower parameter boundary value. The relative sensitivity
percentages are then calculated using equations 26 and 27, where $n_{tot}$ is
the total number of parameters listed in tables 1 and 2.

$SP_{P,i}=100\%\cdot\frac{SI_{P,i}}{\sum_{j=1}^{n_{tot}}SI_{P,j}}$ (26)

$SP_{V,i}=100\%\cdot\frac{SI_{V,i}}{\sum_{j=1}^{n_{tot}}SI_{V,j}}$ (27)

The parameters which describes 99% of the variations in LV volume and pressure
were used in the optimisation algorithm, as will be discussed in section 3.1.

###  2.5 Case studies

To showcase the capabilities of the proposed optimisation model, three case
studies are investigated. The synthetic data generated for case study 1, is
for healthy heart valves using the nominal parameters listed in tables 1 and
2. To recreate clinical data, normally distributed noise [14] was added to the
synthetically generated arterial pressure, chamber volumes and heart valve
flow data which forms the set of measurements used to optimise the model
parameters. The standard deviation used for noise generation was set to 3% of
the mean value for the volume and flow rate measurements and 1% for the
arterial pressure measurement. Case study 2 used the nominal parameters for
data generation but with the aortic valve maximum opening fraction set to
0.222 which corresponds to a valvular flow area of 1 [$cm^{2}$]. Case study 3,
is for a regurgitant mitral valve. For this case, the minimum flow area
opening fraction of the mitral valve was set to 0.05. For each of the case
studies, the three datasets shown in table 3 was used in the optimisation
model to gauge which resulted in the best parameter estimation accuracy.
Figure 3 below shows the synthetically generated chamber volumes, arterial
pressures and valve flow rates generated for the three case studies.

Figure 3: Case studies synthetic data

##  3 Results and discussion

###  3.1 Sensitivity analysis

Figure 4 shows the results of the sensitivity analysis performed using both
automatic differentiation and finite differences to calculate the relevant
gradients for case study 1 (healthy heart). From the autodiff sensitivity
analysis results, it is seen that the LV diastolic elastance parameter has the
largest effect on both LV mean pressure (37.8%) and volume (71.2%). The
parameters with the second and third largest effects on the mean LV pressure
is aortic sinus compliance (18.5%) and LV systolic elastance (12.7%). For the
mean LV volume the second and third rated parameters are LV systolic elastance
(8%) and aortic sinus compliance (5%). Using the sensitivity analysis results,
12 parameters were identified that explained 99% of the variation in mean LV
pressure and volume. The selected 12 optimisation parameters contained in
$\bar{\theta}$ along with the initial systemic and pulmonary venous pressures
and the pulmonary sinus and arterial pressures, can be seen in equation 28.

$\begin{gathered}\bar{\theta}=[E_{lv,s},E_{lv,d},E_{rv,s},E_{rv,d},E_{la,max},E_{la,min},V_{lv,0},C_{AS},\\\
R_{AS},C_{SAT},R_{SAT},R_{PAT},P_{PVN}^{init},P_{SVN}^{init},P_{PS}^{init},P_{PAT}^{init}]\end{gathered}$
(28)

The above sensitivity analysis was also performed for case studies 2 and 3,
and it was found that the algorithm identified the same 12 parameters, albeit,
with a slightly different distribution of sensitivity percentages.

Figure 4: Sensitivity analysis results using autodiff and finite differences

Furthermore, the results in figure 4, show that the parameter sensitivity
percentages calculated for the LV, LA and RV using finite differences are
similar to the values calculated using automatic-differentiation. Conversely,
the finite differences method over- and under-predicts the importance of the
RA parameters and various pulmonary and systemic loop parameters, highlighting
the errors that could be made by approximating gradients using the finite
differences approach.

###  3.2 Case studies

For each of the three case studies discussed in section 2.5, the unknown model
parameters $\bar{\theta}$ are optimised using three optimisation approaches.
The first approach, applies Adam optimisation and forward-mode automatic
differentiation for 1000 training iterations. The second approach uses a
combination of the Adam and LBFGS optimisers, a standard technique used in
scientific machine learning [32], [33]. The model is first trained for 50
iterations using the Adam optimiser after which the second-order optimiser is
utilised. The third approach simply applies the Adam optimiser for 1000
iterations using central finite differences to approximate the required
gradients. For each of the case studies and optimisation approaches the effect
of synthetic dataset contents (available measurements) was also investigated,
by optimising the model parameters using the datasets listed in table 3. For
each of the optimisation runs the initial model parameter vector was kept
constant so as to compare the performance of the different optimisation
permutations.

$\bar{\theta}_{init}=[1,0.01,1,0.01,0.5,0.49,0.1,0.01,1,1,0.1,10,10,25,25]$
(29)

####  3.2.1 Case study 1

Figure 5, shows the mean squared error (MSE) training history for the three
datasets and optimisation approaches using case study 1 model settings and
measurement data. The results indicate that the finite difference optimisation
method (FDM) yields the largest errors compared to the other optimisation
approaches. Note FDM also uses the Adam optimiser but the gradients are
approximated using finite differences. Comparing the training history of the
Adam and hybrid Adam-LBFGS optimisation models (both using forward-mode
automatic differentiation), shows that the latter has superior performance, in
that it results in lower overall MSEs (especially for dataset 2) and converges
in a fewer number of iterations. For the Adam-LBFGS approach, the rapid
reduction of the MSE after the initial 50 iterations can clearly be observed.
It should be noted that for datasets 1 and 3, the Adam optimiser approach
eventually yields MSE values close to the Adam-LBFGS values.

Figure 5: Mean squared error optimisation history for case study 1 and all
datasets

To compare the computational efficiency of the different approaches the total
run times are recorded. The optimisations were performed using a 2020 Apple
MacBook Pro M1. The FDM model calculated 1000 iterations in 230 seconds, the
Adam approach completed the same number of iterations in 71 seconds and the
Adam-LBFGS approach completed its ± 200 iterations in approximately 40
seconds.

To compare the accuracy of the different permutations (optimisation approaches
and datasets), the absolute percentage errors (APEs) between the true nominal
parameters and the optimised values are calculated. Equation 30 shows the
formula used to calculate the APE for the $i^{th}$ parameter, where
$\theta_{i}$ is the optimised parameter and $\theta_{i,true}$ the true nominal
parameter.

$APE=100\%\cdot\frac{\left|\theta_{i}-\theta_{i,true}\right|}{\theta_{i,true}}$
(30)

The absolute and mean absolute percentage error (MAPE) results for case study
1 with the different optimisation approaches and datasets are shown in table
4. Note the shorthand for datasets, D1, D2 and D3 corresponds to datasets 1,2
and 3 respectively. In the tabular results, it is observed that both the Adam
and Adam-LBFGS approaches are capable of estimating the heart parameters
(first seven row entries) with APEs between 0-10%, excluding the unstressed LV
volume, whereas the FDM approach struggles to accurately estimate these
values.

Table 4: Absolute percentage errors per parameter for case study 1 and different synthetic datasets Par. |  D1+ ADAM |  D1+ LBFGS |  D1+ FDM |  D2+ ADAM |  D2+ LBFGS |  D2+ FDM |  D3+ ADAM |  D3+ LBFGS |  D3+ FDM
---|---|---|---|---|---|---|---|---|---
$E_{lv,s}$ |  3.01 |  0.83 |  2.84 |  9.79 |  0.06 |  7.46 |  6.37 |  0.23 |  2.34
$E_{lv,d}$ |  7.11 |  4.67 |  12.60 |  8.25 |  3.05 |  1.37 |  4.22 |  5.17 |  4.26
$E_{rv,s}$ |  6.68 |  5.71 |  15.24 |  1.06 |  1.79 |  3.42 |  0.89 |  2.54 |  31.66
$E_{rv,d}$ |  0.48 |  0.52 |  0.39 |  0.32 |  8.37 |  12.04 |  0.29 |  0.77 |  2.02
$E_{la,max}$ |  8.01 |  5.10 |  9.39 |  8.58 |  5.05 |  19.29 |  5.86 |  6.03 |  5.51
$E_{la,min}$ |  6.65 |  7.29 |  8.13 |  6.30 |  3.41 |  15.11 |  4.70 |  7.56 |  7.68
$V_{lv,0}$ |  18.78 |  33.74 |  61.13 |  5.89 |  39.26 |  23.83 |  35.22 |  38.41 |  63.23
$C_{AS}$ |  54.17 |  27.45 |  81.37 |  178.93 |  2.90 |  71.83 |  1.66 |  4.55 |  87.34
$R_{AS}$ |  77.08 |  42.60 |  188.27 |  113.63 |  3.89 |  57.03 |  96.30 |  6.73 |  44.96
$C_{SAT}$ |  8.59 |  7.91 |  15.60 |  21.02 |  7.16 |  16.69 |  8.93 |  7.54 |  7.95
$R_{SAT}$ |  7.74 |  8.49 |  13.17 |  15.03 |  10.95 |  15.43 |  8.34 |  9.51 |  3.90
$R_{PAT}$ |  24.09 |  1.70 |  29.46 |  6.79 |  4.17 |  14.28 |  12.53 |  3.32 |  33.36
$P_{PVN}^{init}$ |  2.81 |  2.78 |  2.99 |  2.41 |  5.80 |  5.80 |  2.95 |  2.86 |  5.73
$P_{SVN}^{init}$ |  3.62 |  9.28 |  5.50 |  4.50 |  2.86 |  0.64 |  2.03 |  2.81 |  5.84
$P_{PS}^{init}$ |  9.12 |  7.83 |  21.93 |  0.50 |  4.97 |  7.10 |  1.69 |  2.03 |  30.84
MAPE |  15.86 |  11.06 |  31.20 |  25.53 |  6.91 |  18.09 |  12.80 |  6.67 |  22.44

The two configurations with the lowest overall MAPEs are D2 and D3 optimised
using the Adam-LBFGS approach. For both these configurations the $V_{lv,0}$
parameter have error percentages above 30%, but studying the sensitivity
analysis results, it can be seen that the unstressed LV volume has a relative
small effect on the ventricular performance for the analysed range. The result
of this insensitivity to the LV unstressed volume can be seen in figure 6
below, where it is shown that the Adam-LBFGS model configurations using D2 and
D3 can accurately infer the observed arterial pressure and heart chamber
volume waveforms and the unobserved LV pressure waveform by solving the LPM
equations with the optimised parameters. Furthermore, it is seen from table 4,
that although the unstressed volume error is large, the remaining parameter
estimate errors are relatively low. For sake of brevity the other results for
case study 1 are not provided.

Figure 6: Captured heart chamber volumes, arterial pressure and inferred LV
pressure using Adam and Adam-LBFGS models for case study 1 using datasets 2
and 3 (markers – synthetic measurement data, line – simulated using optimised
parameters)

####  3.2.2 Case studies 2 and 3

Seeing as the results generated using D2 and D3 clearly outperform the D1
results for case study 1, the remainder of this section will only use D2 and
D3 for parameter estimations. Figure 7 shows the training history for case
studies 2 (aortic stenosis) and 3 (mitral regurgitation) using D2 and D3 along
with the mentioned optimisation approaches. The results again show that the
Adam-LBFGS approach outperforms the other two optimisers in terms of converged
MSE values. The converged MSE values generated using D2 for both case studies
are approximately $10^{-3}$, whereas, the MSE values achieved using D3 for
model training are approximately $10^{-2}$. This leads to the conclusion that
the addition of the heart chamber volume waveforms in the measurement data
increases the resultant MSEs. This does not necessarily mean that using D2
will produce more accurate parameter estimates than using D3 as will be
discussed in the next paragraph.

Figure 7: Mean squared optimisation error history for case studies 2 and 3
using datasets 2 and 3

Table 5 contains the MAPEs between the true and estimated model parameters for
case studies 2 and 3. The results show that the Adam-LBFGS optimiser approach
for both cases produces smaller errors when compared to the Adam-only
optimised parameter values. For case 2 utilising the Adam-LBFGS approach, the
results achieved using only the arterial pressure and valvular flow data (D2)
yielded smaller errors when compared to the results when using D3. Conversely,
for case 3, the addition of heart chamber volume data (D3) resulted in the
smallest MAPE. Therefore, using D2 for the aortic stenosis case and D3 for the
mitral regurgitation case along with the Adam-LBFGS optimiser resulted in the
best parameter estimates.

Table 5: Mean absolute percentage errors for case studies 2 and 3 using datasets 2 and 3 | Case study 2 | Case study 3
---|---|---
|  Adam |  Adam + LBFGS |  FDM |  Adam |  Adam + LBFGS |  FDM
Dataset 2 |  17.50 |  7.68 |  17.09 |  8.01 |  6.79 |  88.14
Dataset 3 |  18.95 |  14.14 |  21.03 |  7.05 |  6.16 |  76.17

To further investigate the effects of measurement data and valvular diseases
on the heart parameter estimates, bar plots of the individual APEs are plotted
in figure 8 for the different optimisation approaches. The results show that
the FDM struggles to accurately estimate the model parameters for case study
3, and performs similar to the Adam optimised model for case study 2 which
aligns with training history observed in figure 7. Studying the Adam and Adam-
LBFGS error plots for both case studies, it is seen that using D2 for model
training, results in smaller LV elastance parameter errors when compared to
the error results generated when using D3, but generally using D3 for model
training yields smaller errors for the remaining parameters shown in figure 8
along with the various vasculature parameters (not shown). Therefore, we can
conclude that the addition of the volume waveform information into the
measurement dataset does not produce more accurate LV elastance parameter
estimates but does increase the accuracy of the remaining parameter estimates
at the cost of a slight increase in LV elastance parameter errors.

Figure 8: Heart parameter absolute percentage errors for case studies 2 and 3
(X-axis legend see figure 4)

To visualise the effect of using the different measurement datasets, the
predicted and actual volume, pressure and flow rate wave forms are plotted in
figures 9 and 10 for parameter sets estimated using D2 and D3 with the Adam-
LBFGS optimiser. It is seen in figure 9 that the predicted parameters can
clearly replicate the heart chamber volume waveforms, along with the arterial
pressure and valvular flow rates for both case studies. For the unobserved LV
pressure it is seen that the model slightly under predicts the pressure for
case 2 and slightly overpredicts the LV pressure for case 3.

Figure 9: Predicted heart chamber volumes, arterial pressure and valvular flow
rates using dataset 2 (markers – synthetic measurement data, line – simulated
using optimised parameters)

In figure 10, similar to the results in figure 9, it is seen that the LPM with
the optimised parameters using D3 can with relative accuracy replicate the
actual waveforms. It is seen that the predicted LV pressure accurately tracks
the actual LV pressure, and in contrast to the results in figure 9, does not
over- or under-shoot the true waveform. This can be ascribed to the fact that
the addition of heart chamber volume waveforms in the measurement data, leads
to the optimiser more accurately capturing the various compartment
resistances, capacitances and blood inertias which make up the ventricular
afterload. Therefore, if the equipment is available to accurately measure the
atrial and ventricular volume changes it is recommended to use the volume data
in addition to the flow rate and pressure data when optimising the model
parameters to accurately capture the ventricular afterload. That being said,
using only the flow rate and pressure data, based on the current findings,
would still yield accurate elastance parameter estimates.

Figure 10: Predicted heart chamber volumes, arterial pressure and valvular
flow rates using dataset 3 (markers – synthetic measurement data, line –
simulated using optimised parameters)

##  4 Limitations

In the present work, it is assumed that the lumped parameter/0D modelling
approach is sufficient to capture the dynamics of an actual cardiovascular
system. Furthermore, the heart valve flow coefficients were kept at constant
values taken from literature, but in reality these flow coefficients would
change based on the varying flow Reynolds number in the valve. Another
limitation of the current work, is that the sensitivity analysis results are
dependent on the selected parameter ranges, which in the current work is not
based on values found in literature and rather selected to cover a wide range
of vasculature impedances and heart chamber elastances.

##  5 Conclusions

The goal of the present work was to demonstrate the ability of computational
parameter estimation algorithms to estimate LV elastances for diseased heart
valve cases. It was shown using synthetically generated measurement data, that
the use of automatic differentiation over finite differences to calculate the
required gradients for cardiovascular parameter estimation, significantly
improves the overall accuracy of the optimised parameters as well as the
computational efficiency. Various gradient-based optimisers were also
evaluated and it was found that a combination of Adam and LBFGS optimisers
produced the best results for the analysed configurations. Using the hybrid
Adam-LBFGS optimiser, important parameters related to LV performance was
estimated for cases with aortic stenosis and mitral regurgitation using
different collections of measurements. It was found that using arterial
pressure and valvular flow rate waveforms along with the proposed optimiser,
the approach was capable of predicting the LV elastances within approximately
2% of the actual values. Adding the ventricular and atrial volume waveforms to
the training data increased the LV elastance value errors to approximately 5%,
but the additional information meant the model was capable of more accurately
capturing the various vasculature parameters. The combination of slightly
higher LV elastance and lower vasculature parameter errors, resulted in
simulation results that more closely track the actual LV pressure waveforms.
Future work entails validating the proposed approach with actual patient data
with known LV elastances.

## References

  * [1] R. BeLue, T. A. Okoror, J. Iwelunmor, K. D. Taylor, A. N. Degboe, C. Agyemang, G. Ogedegbe, An overview of cardiovascular risk factor burden in sub-Saharan African countries: A socio-cultural perspective, Globalization and Health 5 (2009) 1–12.  doi:10.1186/1744-8603-5-10.
  * [2] A. K. Keates, A. O. Mocumbi, M. Ntsekhe, K. Sliwa, S. Stewart, Cardiovascular disease in Africa: Epidemiological profile and challenges, Nature Reviews Cardiology 14 (5) (2017) 273–293.  doi:10.1038/nrcardio.2017.19.
  * [3] M. F. Yuyun, K. Sliwa, A. P. Kengne, A. O. Mocumbi, G. Bukhman, Cardiovascular diseases in sub-saharan Africa compared to high-income countries: An epidemiological perspective, Global Heart 15 (1) (2020) 1–18.  doi:10.5334/GH.403.
  * [4] C. Serbia, S. B. Germany, J. B. Germany, L. C. Germany, N. Freemantle, U. Kingdom, M. G. France, P. Ju, C. T. France, W. W. Poland, E. S. C. Eacts, 2021 ESC / EACTS Guidelines for the management of valvular heart disease Developed by the Task Force for the management of valvular heart disease of the European Society of Cardiology ( ESC ) and the, European Heart Journal (2021) 1–72.
  * [5] E. Gayat, V. Mor-Avi, L. Weinert, C. Yodwut, R. M. Lang, Noninvasive quantification of left ventricular elastance and ventricular-arterial coupling using three-dimensional echocardiography and arterial tonometry, American Journal of Physiology - Heart and Circulatory Physiology 301 (5) (2011) 1916–1923.  doi:10.1152/ajpheart.00760.2011.
  * [6] D. A. Kass, T. Yamazaki, D. Burkhoff, W. L. Maughan, K. Sagawa, Determination of left ventricular end-systolic pressure-volume relationships by the conductance (volume) catheter technique, Circulation 73 (3) (1986) 586–595.  doi:10.1161/01.CIR.73.3.586.
  * [7] C. H. Chen, B. Fetics, E. Nevo, C. E. Rochitte, K. R. Chiou, P. Y. A. Ding, M. Kawaguchi, D. A. Kass, Noninvasive single-beat determination of left ventricular end-systolic elastance in humans, Journal of the American College of Cardiology 38 (7) (2001) 2028–2034.  doi:10.1016/S0735-1097(01)01651-5.
  * [8] V. Bikia, D. Adamopoulos, S. Pagoulatou, G. Rovas, N. Stergiopulos, AI-Based Estimation of End-Systolic Elastance From Arm-Pressure and Systolic Time Intervals, Frontiers in Artificial Intelligence 4 (April) (2021) 1–11.  doi:10.3389/frai.2021.579541.
  * [9] F. Huang, S. Ying, On-line parameter identification of the lumped arterial system model: A simulation study, PLoS ONE 15 (7 July) (2020) 1–17.  doi:10.1371/journal.pone.0236012.
URL http://dx.doi.org/10.1371/journal.pone.0236012

  * [10] Mathworks, fmincon (2022).
URL https://www.mathworks.com/help/optim/ug/fmincon.html

  * [11] X. Zhang, D. Wu, F. Miao, H. Liu, Y. Li, Personalized Hemodynamic Modeling of the Human Cardiovascular System: A Reduced-Order Computing Model, IEEE Transactions on Biomedical Engineering 67 (10) (2020) 2754–2764.  doi:10.1109/TBME.2020.2970244.
  * [12] A. L. Colunga, K. G. Kim, N. P. Woodall, T. F. Dardas, J. H. Gennari, M. S. Olufsen, B. E. Carlson, Deep phenotyping of cardiac function in heart transplant patients using cardiovascular system models, Journal of Physiology 598 (15) (2020) 3203–3222.  arXiv:1812.11857, doi:10.1113/JP279393.
  * [13] Z. Keshavarz-Motamed, A diagnostic, monitoring, and predictive tool for patients with complex valvular, vascular and ventricular diseases, Scientific Reports 10 (1) (2020) 1–19.  doi:10.1038/s41598-020-63728-8.
URL http://dx.doi.org/10.1038/s41598-020-63728-8

  * [14] N. L. Bjordalsbakke, J. T. Sturdy, D. R. Hose, L. R. Hellevik, Parameter estimation for closed-loop lumped parameter models of the systemic circulation using synthetic data, Mathematical Biosciences 343 (October 2021) (2022) 108731.  doi:10.1016/j.mbs.2021.108731.
URL https://doi.org/10.1016/j.mbs.2021.108731

  * [15] M. A. Branch, T. Coleman, Y. Li, A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems, SIAM Journal on Scientific Computing 21 (1) (1999) 1–23.
  * [16] SciPy Organisation, SciPy 1.8.0: User guide (2021).
URL
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html

  * [17] A. L. Marsden, Optimization in cardiovascular modeling, Annual Review of Fluid Mechanics 46 (September) (2014) 519–546.  doi:10.1146/annurev-fluid-010313-141341.
  * [18] M. J. Kochenderfer, T. Wheeler, Algorithms for Optimization, 1st Edition, MIT Press, Cambridge, 2019.
  * [19] C. Rackauckas, Q. Nie, DifferentialEquations.jl, Journal of Open Research Software 5 (1) (2017) 15.
URL https://diffeq.sciml.ai/stable/

  * [20] J. Revels, M. Lubin, T. Papamarkou, Forward-Mode Automatic Differentiation in Julia, arXiv:1607.07892.
  * [21] P. Morgensen, A. Riseth, Optim: A mathematical optimization package for Julia, Journal of Open Source Software 3 (24) (2018) 615.
  * [22] M. Innes, Flux: Elegant Machine Learning with Julia, Journal of Open Source Software.
  * [23] H. Suga, K. Sagawa, A. A. Shoukas, Load independence of the instantaneous pressure-volume ratio of the canine left ventricle and effects of epinephrine and heart rate on the ratio., Circulation research 32 (3) (1973) 314–322.  doi:10.1161/01.RES.32.3.314.
  * [24] T. Korakianitis, Y. Shi, A concentrated parameter model for the human cardiovascular system including heart valve dynamics and atrioventricular interaction, Medical Engineering and Physics 28 (7) (2006) 613–628.  doi:10.1016/j.medengphy.2005.10.004.
  * [25] S. Bozkurt, Mathematical modeling of cardiac function to evaluate clinical cases in adults and children, PLoS ONE 14 (10) (2019) 1–20.  doi:10.1371/journal.pone.0224663.
URL http://dx.doi.org/10.1371/journal.pone.0224663

  * [26] T. Korakianitis, Y. Shi, Numerical simulation of cardiovascular dynamics with healthy and diseased heart valves, Journal of Biomechanics 39 (11) (2006) 1964–1982.  doi:10.1016/j.jbiomech.2005.06.016.
  * [27] H. Tang, Z. Dai, M. Wang, B. Guo, S. Wang, J. Wen, T. Li, Lumped-Parameter Circuit Platform for Simulating Typical Cases of Pulmonary Hypertensions from Point of Hemodynamics, Journal of Cardiovascular Translational Research 13 (5) (2020) 826–852.  doi:10.1007/s12265-020-09953-y.
  * [28] C. Tsitouras, Runge-Kutta pairs of order 5(4) satisfying only the first column simplifying assumption, Computers and Mathematics with Applications 62 (2) (2011) 770–775.  doi:10.1016/j.camwa.2011.06.002.
URL http://dx.doi.org/10.1016/j.camwa.2011.06.002

  * [29] R. Laubscher, J. Liebenberg, P. Herbst, Dynamic simulation of aortic valve stenosis using a lumped parameter cardiovascular system model with flow regime dependent valve pressure loss characteristicsarXiv:2204.03701.
URL http://arxiv.org/abs/2204.03701

  * [30] D. P. Kingma, J. L. Ba, Adam: A method for stochastic optimization, arXiv preprint (2015) 1–15arXiv:arXiv:1412.6980v9.
  * [31] J. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Mathematics of Computation 35 (151) (1980) 773–782.
  * [32] R. Laubscher, P. Rousseau, Application of a mixed variable physics-informed neural network to solve the incompressible steady-state and transient mass, momentum, and energy conservation equations for flow over in-line heated tubes, Applied Soft Computing (xxxx) (2021) 108050\.  doi:10.1016/j.asoc.2021.108050.
URL https://doi.org/10.1016/j.asoc.2021.108050

  * [33] S. Markidis, The Old and the New: Can Physics-Informed Deep-Learning Replace Traditional Linear Solvers? 1 (1) (2021) 1–20.  arXiv:2103.09655, doi:10.3389/fdata.2021.669097.
URL http://arxiv.org/abs/2103.09655


pursuit of understanding fundamental limits of light-matter interaction due to
their universal applicability.

###### Acknowledgements.

The research leading to the presented results has received funding from the
European Research Council under the European Union’s Seventh Framework ERC
Grant Agreement No. 615613. We also express gratitude to the DFG (German
Science Foundation) for continuous funding of our research via the projects
RE2974/5-1, RE2974/9-1, RE2974/12-1, and CRC 787. AC gratefully acknowledges
support from the Deutsche Forschungsgemeinschaft (DFG) through the project B1
of the SFB 910 and from the European Unions Horizon 2020 research and
innovation program under the SONAR grant Agreement No. 734690 and Andreas
Knorr for fruitful discussions.

## References

  * Michler  _et al._ (2000) P. Michler, A. Kiraz, C. Becher, W. Schoenfeld, P. Petroff, L. Zhang, E. Hu, and A. Imamoglu, science 290, 2282 (2000).
  * Yuan  _et al._ (2002) Z. Yuan, B. E. Kardynal, R. M. Stevenson, A. J. Shields, C. J. Lobo, K. Cooper, N. S. Beattie, D. A. Ritchie, and M. Pepper, science 295, 102 (2002).
  * Akopian  _et al._ (2006) N. Akopian, N. Lindner, E. Poem, Y. Berlatzky, J. Avron, D. Gershoni, B. Gerardot, and P. Petroff, Physical Review Letters 96, 130501 (2006).
  * Hafenbrak  _et al._ (2007) R. Hafenbrak, S. Ulrich, P. Michler, L. Wang, A. Rastelli, and O. Schmidt, New Journal of Physics 9, 315 (2007).
  * Santori  _et al._ (2002) C. Santori, D. Fattal, J. Vučković, G. S. Solomon, and Y. Yamamoto, Nature 419, 594 (2002).
  * Reithmaier  _et al._ (2004) J. P. Reithmaier, G. Sek, A. Löffler, C. Hofmann, S. Kuhn, S. Reitzenstein, L. Keldysh, V. Kulakovskii, T. Reinecke, and A. Forchel, Nature 432, 197 (2004).
  * Woggon (1997) U. Woggon,  _Optical properties of semiconductor quantum dots_ (Springer, 1997).
  * Bimberg  _et al._ (1999) D. Bimberg, M. Grundmann, and N. N. Ledentsov,  _Quantum dot heterostructures_ (John Wiley & Sons, 1999).
  * Kiraz  _et al._ (2004) A. Kiraz, M. Atatüre, and A. Imamoğlu, Physical Review A 69, 032305 (2004).
  * Michler (2009) P. Michler,  _Single semiconductor quantum dots_ , Vol. 231 (Springer, 2009).
  * Michler (2017) P. Michler,  _Quantum Dots for Quantum Information Technologies_ , Vol. 237 (Springer, 2017).
  * Jahnke (2012) F. Jahnke,  _Quantum optics with semiconductor nanostructures_ (Elsevier, 2012).
  * Shields (2010) A. J. Shields, in  _Nanoscience And Technology: A Collection of Reviews from Nature Journals_ (World Scientific, 2010) pp. 221–229.
  * Yu-Ming  _et al._ (2013) H. Yu-Ming, H. Yu, W. Yu-Jia, W. Dian, A. Mete, S. Christian, H. Sven, K. Martin, L. Chao-Yang, and P. Jian-Wei, Nature Nanotechnology 8, 213 (2013).
  * Peter  _et al._ (2004) E. Peter, J. Hours, P. Senellart, A. Vasanelli, A. Cavanna, J. Bloch, and J. M. Gérard, Phys. Rev. B 69, 041307 (2004).
  * Mendach  _et al._ (2006) S. Mendach, R. Songmuang, S. Kiravittaya, A. Rastelli, M. Benyoucef, and O. G. Schmidt, Applied Physics Letters 88, 111120 (2006).
  * Schmidt (2007) O. G. Schmidt,  _Lateral alignment of epitaxial quantum dots_ (Springer Science & Business Media, 2007).
  * Zhang  _et al._ (2015) J. Zhang, J. S. Wildmann, F. Ding, R. Trotta, Y. Huo, E. Zallo, D. Huber, A. Rastelli, and O. G. Schmidt, Nature communications 6, 10067 (2015).
  * Schlehahn  _et al._ (2018) A. Schlehahn, S. Fischbach, R. Schmidt, A. Kaganskiy, A. Strittmatter, S. Rodt, T. Heindel, and S. Reitzenstein, Scientific reports 8, 1340 (2018).
  * Yoshie  _et al._ (2004) T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova, H. Gibbs, G. Rupper, C. Ell, O. Shchekin, and D. Deppe, Nature 432, 200 (2004).
  * Bose  _et al._ (2014) R. Bose, T. Cai, K. R. Choudhury, G. S. Solomon, and E. Waks, Nature Photonics 8, 858 (2014).
  * Senellart  _et al._ (2017) P. Senellart, G. Solomon, and A. White, Nature nanotechnology 12, 1026 (2017).
  * Förstner  _et al._ (2003) J. Förstner, C. Weber, J. Danckwerts, and A. Knorr, physica status solidi (b) 238, 419 (2003).
  * Krummheuer  _et al._ (2005) B. Krummheuer, V. M. Axt, T. Kuhn, I. D’Amico, and F. Rossi, Phys. Rev. B 71, 235329 (2005).
  * Quilter  _et al._ (2015) J. H. Quilter, A. J. Brash, F. Liu, M. Glässl, A. M. Barth, V. M. Axt, A. J. Ramsay, M. S. Skolnick, and A. M. Fox, Phys. Rev. Lett. 114, 137401 (2015).
  * Müller  _et al._ (2014) M. Müller, S. Bounouar, K. D. Jöns, M. Glässl, and P. Michler, Nature Photonics 8, 224 (2014).
  * Reiter  _et al._ (2014) D. Reiter, T. Kuhn, M. Glässl, and V. Axt, Journal of Physics: Condensed Matter 26, 423203 (2014).
  * Carmele  _et al._ (2010) A. Carmele, M. Richter, W. W. Chow, and A. Knorr, Phys. Rev. Lett. 104, 156801 (2010).
  * Carmele  _et al._ (2013a) A. Carmele, A. Knorr, and F. Milde, New Journal of Physics 15, 105024 (2013a).
  * Besombes  _et al._ (2001a) L. Besombes, K. Kheng, L. Marsal, and H. Mariette, Physical Review B 63, 155307 (2001a).
  * Borri  _et al._ (2001) P. Borri, W. Langbein, S. Schneider, U. Woggon, R. L. Sellin, D. Ouyang, and D. Bimberg, Physical Review Letters 87, 157401 (2001).
  * Galland  _et al._ (2008) C. Galland, A. Högele, H. E. Türeci, and A. Imamoğlu, Physical review letters 101, 067402 (2008).
  * Vagov  _et al._ (2007) A. Vagov, M. Croitoru, V. M. Axt, T. Kuhn, and F. Peeters, Physical review letters 98, 227403 (2007).
  * Madsen  _et al._ (2011) K. H. Madsen, S. Ates, T. Lund-Hansen, A. Löffler, S. Reitzenstein, A. Forchel, and P. Lodahl, Physical review letters 106, 233601 (2011).
  * Carmichael (2009) H. J. Carmichael,  _Statistical Methods in Quantum Optics 1_ (Springer Science & Business Media, 2009).
  * Scully and Zubairy (1999) M. O. Scully and M. S. Zubairy, “Quantum optics,” (1999).
  * Breuer  _et al._ (2002) H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Gardiner  _et al._ (2004) C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ , Vol. 56 (Springer Science & Business Media, 2004).
  * Breuer  _et al._ (2016) H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016).
  * Cerrillo and Cao (2014) J. Cerrillo and J. Cao, Phys. Rev. Lett. 112, 110401 (2014).
  * de Vega and Alonso (2017) I. de Vega and D. Alonso, Rev. Mod. Phys. 89, 015001 (2017).
  * Michler (2003) P. Michler,  _Single quantum dots: Fundamentals, applications and new concepts_ , Vol. 90 (Springer Science & Business Media, 2003).
  * Stock  _et al._ (2011) E. Stock, M.-R. Dachner, T. Warming, A. Schliwa, A. Lochmann, A. Hoffmann, A. I. Toropov, A. K. Bakarov, I. A. Derebezov, M. Richter, V. A. Haisler, A. Knorr, and D. Bimberg, Phys. Rev. B 83, 041304 (2011).
  * Strauß  _et al._ (2018) M. Strauß, M. Hohn, C. Schneider, S. Höfling, J. Wolters, and S. Reitzenstein, arXiv preprint arXiv:1805.06357 (2018).
  * Ramsay  _et al._ (2010a) A. Ramsay, T. Godden, S. Boyle, E. M. Gauger, A. Nazir, B. W. Lovett, A. Fox, and M. Skolnick, Physical Review Letters 105, 177402 (2010a).
  * Kaldewey  _et al._ (2017) T. Kaldewey, S. Lüker, A. V. Kuhlmann, S. R. Valentin, J.-M. Chauveau, A. Ludwig, A. D. Wieck, D. E. Reiter, T. Kuhn, and R. J. Warburton, Physical Review B 95, 241306 (2017).
  * Bounouar  _et al._ (2015) S. Bounouar, M. Müller, A. Barth, M. Glässl, V. Axt, and P. Michler, Physical Review B 91, 161302 (2015).
  * Ardelt  _et al._ (2014) P.-L. Ardelt, L. Hanschke, K. A. Fischer, K. Müller, A. Kleinkauf, M. Koller, A. Bechtold, T. Simmet, J. Wierzbowski, H. Riedl, _et al._ , Physical Review B 90, 241404 (2014).
  * Iles-Smith and Nazir (2016) J. Iles-Smith and A. Nazir, Optica 3, 207 (2016).
  * Weiler  _et al._ (2012) S. Weiler, A. Ulhaq, S. M. Ulrich, D. Richter, M. Jetter, P. Michler, C. Roy, and S. Hughes, Phys. Rev. B 86, 241304 (2012).
  * Hohenester (2010) U. Hohenester, Phys. Rev. B 81, 155303 (2010).
  * Monniello  _et al._ (2013) L. Monniello, C. Tonin, R. Hostein, A. Lemaitre, A. Martinez, V. Voliotis, and R. Grousson, Physical Review Letters 111, 026403 (2013).
  * Ulrich  _et al._ (2011) S. Ulrich, S. Ates, S. Reitzenstein, A. Löffler, A. Forchel, and P. Michler, Physical review letters 106, 247402 (2011).
  * Roy and Hughes (2011) C. Roy and S. Hughes, Physical review letters 106, 247403 (2011).
  * Gardiner and Zoller (1991) C. Gardiner and P. Zoller,  _Quantum Noise_ (Springer, Berlin Heidelberg New York, 1991).
  * Koch  _et al._ (1993) S. W. Koch _et al._ ,  _Semiconductor quantum dots_ , Vol. 2 (World Scientific, 1993).
  * Grundmann  _et al._ (1995) M. Grundmann, O. Stier, and D. Bimberg, Phys. Rev. B 52, 11969 (1995).
  * Stier  _et al._ (1999) O. Stier, M. Grundmann, and D. Bimberg, Phys. Rev. B 59, 5688 (1999).
  * Schliwa  _et al._ (2007) A. Schliwa, M. Winkelnkemper, and D. Bimberg, Phys. Rev. B 76, 205324 (2007).
  * Chow and Jahnke (2013) W. W. Chow and F. Jahnke, Progress in quantum electronics 37, 109 (2013).
  * Haug and Koch (2009) H. Haug and S. W. Koch,  _Quantum Theory of the Optical and Electronic Properties of Semiconductors: Fivth Edition_ (World Scientific Publishing Company, 2009).
  * Brandes (2005) T. Brandes, physics reports 408, 315 (2005).
  * Sakaki (1992) H. Sakaki, Surface Science 267, 623 (1992).
  * Schwab  _et al._ (2000) K. Schwab, E. Henriksen, J. Worlock, and M. L. Roukes, Nature 404, 974 (2000).
  * Gornyi  _et al._ (2005) I. Gornyi, A. Mirlin, and D. Polyakov, Physical review letters 95, 206603 (2005).
  * Sun  _et al._ (2000) C.-K. Sun, J.-C. Liang, and X.-Y. Yu, Physical review letters 84, 179 (2000).
  * Lugli and Goodnick (1987) P. Lugli and S. Goodnick, Physical review letters 59, 716 (1987).
  * Kim  _et al._ (1992) D.-S. Kim, J. Shah, J. Cunningham, T. Damen, W. Schäfer, M. Hartmann, and S. Schmitt-Rink, Physical review letters 68, 1006 (1992).
  * Gammon  _et al._ (1995) D. Gammon, S. Rudin, T. Reinecke, D. Katzer, and C. Kyono, Physical Review B 51, 16785 (1995).
  * Richter  _et al._ (2018) M. Richter, R. Singh, M. Siemens, and S. T. Cundiff, Science advances 4, eaar7697 (2018).
  * Singh  _et al._ (2017) R. Singh, M. Richter, G. Moody, M. E. Siemens, H. Li, and S. T. Cundiff, Phys. Rev. B 95, 235307 (2017).
  * Lim  _et al._ (2017) J. Lim, D. J. Ing, J. Rosskopf, J. Jeske, J. H. Cole, S. F. Huelga, and M. B. Plenio, The Journal of chemical physics 146, 024109 (2017).
  * Wigger  _et al._ (2018) D. Wigger, C. Schneider, S. Gerhardt, M. Kamp, S. Höfling, T. Kuhn, and J. Kasprzak, Optica 5, 1442 (2018).
  * Cassette  _et al._ (2015) E. Cassette, R. D. Pensack, B. Mahler, and G. D. Scholes, Nature communications 6, 6086 (2015).
  * Liu  _et al._ (2018a) A. Liu, D. B. Almeida, W. K. Bae, L. A. Padilha, and S. T. Cundiff, arXiv preprint arXiv:1806.06112 (2018a).
  * Suzuki  _et al._ (2018) T. Suzuki, R. Singh, G. Moody, M. Aßmann, M. Bayer, A. Ludwig, A. D. Wieck, and S. T. Cundiff, Phys. Rev. B 98, 195304 (2018).
  * Christiansen  _et al._ (2017) D. Christiansen, M. Selig, G. Berghäuser, R. Schmidt, I. Niehues, R. Schneider, A. Arora, S. M. de Vasconcellos, R. Bratschitsch, E. Malic, and A. Knorr, Phys. Rev. Lett. 119, 187402 (2017).
  * Malic and Knorr (2013) E. Malic and A. Knorr,  _Graphene and Carbon Nanotubes: Ultrafast Optics and Relaxation Dynamics_ (John Wiley & Sons, 2013).
  * Kira and Koch (2011) M. Kira and S. W. Koch,  _Semiconductor quantum optics_ (Cambridge University Press, 2011).
  * Hohenester  _et al._ (2009) U. Hohenester, A. Laucht, M. Kaniber, N. Hauke, A. Neumann, A. Mohtashami, M. Seliger, M. Bichler, and J. J. Finley, Physical Review B 80, 201311 (2009).
  * Briegel  _et al._ (1998) H.-J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 81, 5932 (1998).
  * H.J. Kimble (2008) H.J. Kimble, Nature 453, 1023 (2008).
  * Mukamel (1995) S. Mukamel,  _Principles of nonlinear optical spectroscopy_ , Vol. 29 (Oxford university press New York, 1995).
  * Axt and Mukamel (1998) V. M. Axt and S. Mukamel, Rev. Mod. Phys. 70, 145 (1998).
  * R. J. Warburton (2013) R. J. Warburton, Nature Materials 12, 483 (2013).
  * Matthiesen  _et al._ (2012) C. Matthiesen, A. N. Vamivakas, and M. Atatüre, Phys. Rev. Lett. 108, 093602 (2012).
  * May and Kühn (2008) V. May and O. Kühn,  _Charge and energy transfer dynamics in molecular systems_ (John Wiley & Sons, 2008).
  * Bourgain  _et al._ (2013) R. Bourgain, J. Pellegrino, S. Jennewein, Y. R. Sortais, and A. Browaeys, Optics letters 38, 1963 (2013).
  * Sondermann and Leuchs (2013) M. Sondermann and G. Leuchs, Journal of the European Optical Society - Rapid publications 8 (2013).
  * Palatchi  _et al._ (2014) C. Palatchi, J. M. Dahlstr m, A. S. Kheifets, I. A. Ivanov, D. M. Canaday, P. Agostini, and L. F. DiMauro, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 245003 (2014).
  * Kira and Koch (2006) M. Kira and S. Koch, Progress in quantum electronics 30, 155 (2006).
  * Čížek (1969) J. Čížek, Advances in chemical physics , 35 (1969).
  * Krummheuer  _et al._ (2002) B. Krummheuer, V. M. Axt, and T. Kuhn, Phys. Rev. B 65, 195313 (2002).
  * Richter  _et al._ (2009) M. Richter, A. Carmele, A. Sitek, and A. Knorr, Phys. Rev. Lett. 103, 087407 (2009).
  * Rozbicki and Machnikowski (2008) E. Rozbicki and P. Machnikowski, Phys. Rev. Lett. 100, 027401 (2008).
  * Lüker  _et al._ (2017) S. Lüker, T. Kuhn, and D. E. Reiter, Phys. Rev. B 96, 245306 (2017).
  * Wilson-Rae and Imamoğlu (2002) I. Wilson-Rae and A. Imamoğlu, Phys. Rev. B 65, 235311 (2002).
  * Laussy  _et al._ (2012) F. Laussy, E. D. Valle, T. M nchen, A. Laucht, A. Gonzalez-Tudela, M. Kaniber, J. Finley, and C. Tejedor, in  _Quantum Optics with Semiconductor Nanostructures_, Woodhead Publishing Series in Electronic and Optical Materials, edited by F. Jahnke (Woodhead Publishing, 2012) pp. 293 – 331.
  * Laussy  _et al._ (2011) F. P. Laussy, A. Laucht, E. del Valle, J. J. Finley, and J. M. Villas-Bôas, Phys. Rev. B 84, 195313 (2011).
  * del Valle (2010) E. del Valle, Saarbrcken: VDM Verlag (2010).
  * del Valle and Laussy (2010) E. del Valle and F. P. Laussy, Phys. Rev. Lett. 105, 233601 (2010).
  * Kreinberg  _et al._ (2018) S. Kreinberg, T. Grbešić, M. Strauß, A. Carmele, M. Emmerling, C. Schneider, S. Höfling, X. Porte, and S. Reitzenstein, Light, Science & Applications 7 (2018).
  * McCutcheon  _et al._ (2011) D. P. S. McCutcheon, N. S. Dattani, E. M. Gauger, B. W. Lovett, and A. Nazir, Phys. Rev. B 84, 081305 (2011).
  * Manson  _et al._ (2016) R. Manson, K. Roy-Choudhury, and S. Hughes, Phys. Rev. B 93, 155423 (2016).
  * Kabuss  _et al._ (2011a) J. Kabuss, A. Carmele, M. Richter, W. W. Chow, and A. Knorr, physica status solidi (b) 248, 872 (2011a).
  * Axt  _et al._ (1999) V. Axt, M. Herbst, and T. Kuhn, Superlattices and Microstructures 26, 117 (1999).
  * Kabuss  _et al._ (2011b) J. Kabuss, A. Carmele, M. Richter, and A. Knorr, Phys. Rev. B 84, 125324 (2011b).
  * Glässl  _et al._ (2011a) M. Glässl, A. Vagov, S. Lüker, D. E. Reiter, M. D. Croitoru, P. Machnikowski, V. M. Axt, and T. Kuhn, Phys. Rev. B 84, 195311 (2011a).
  * Schollwöck (2011) U. Schollwöck, Annals of Physics 326, 96 (2011).
  * Caldeira and Leggett (1983) A. O. Caldeira and A. J. Leggett, Physica A: Statistical mechanics and its Applications 121, 587 (1983).
  * Vidal (2004) G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * del Pino  _et al._ (2018) J. del Pino, F. A. Y. N. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, Phys. Rev. Lett. 121, 227401 (2018).
  * Strathearn  _et al._ (2018) A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Nature communications 9, 3322 (2018).
  * Droenner  _et al._ (2018) L. Droenner, N. L. Naumann, A. Knorr, and A. Carmele, arXiv preprint arXiv:1801.03342 (2018).
  * Makri (1998) N. Makri, The Journal of Physical Chemistry A 102, 4414 (1998).
  * Cygorek  _et al._ (2017) M. Cygorek, A. M. Barth, F. Ungar, A. Vagov, and V. M. Axt, Physical Review B 96, 201201 (2017).
  * Hopfmann  _et al._ (2015) C. Hopfmann, A. Musiał, M. Strauß, A. M. Barth, M. Glässl, A. Vagov, M. Strauß, C. Schneider, S. Höfling, M. Kamp, V. M. Axt, and S. Reitzenstein, Phys. Rev. B 92, 245403 (2015).
  * Jackson (2012) J. D. Jackson,  _Classical electrodynamics_ (John Wiley & Sons, 2012).
  * Allen and Eberly (1987) L. Allen and J. H. Eberly,  _Optical resonance and two-level atoms_ , Vol. 28 (Courier Corporation, 1987).
  * Gardiner (2009) C. Gardiner,  _Stochastic methods_ , Vol. 4 (springer Berlin, 2009).
  * Besombes  _et al._ (2001b) L. Besombes, K. Kheng, L. Marsal, and H. Mariette, Physical Review B 63, 155307 (2001b).
  * Jakubczyk  _et al._ (2016) T. Jakubczyk, V. Delmonte, S. Fischbach, D. Wigger, D. E. Reiter, Q. Mermillod, P. Schnauber, A. Kaganskiy, J.-H. Schulze, A. Strittmatter, _et al._ , ACS photonics 3, 2461 (2016).
  * Fras  _et al._ (2016) F. Fras, Q. Mermillod, G. Nogues, C. Hoarau, C. Schneider, M. Kamp, S. Höfling, W. Langbein, and J. Kasprzak, Nature Photonics 10, 155 (2016).
  * Madelung (2012) O. Madelung,  _Introduction to solid-state theory_ , Vol. 2 (Springer Science & Business Media, 2012).
  * Pullerits  _et al._ (1995) T. Pullerits, R. Monshouwer, F. van Mourik, and R. van Grondelle, Chemical physics 194, 395 (1995).
  * Cohen-Tannoudji and Guéry-Odelin (2011) C. Cohen-Tannoudji and D. Guéry-Odelin,  _Advances in atomic physics: an overview_ (World Scientific, 2011).
  * Gardiner and Zoller (2015) C. Gardiner and P. Zoller, in  _The Quantum World of Ultra-Cold Atoms and Light Book II: The Physics of Quantum-Optical Devices_ (World Scientific, 2015) pp. 1–524.
  * McCutcheon and Nazir (2010) D. P. McCutcheon and A. Nazir, New Journal of Physics 12, 113042 (2010).
  * Kabuss  _et al._ (2010) J. Kabuss, S. Werner, A. Hoffmann, P. Hildebrandt, A. Knorr, and M. Richter, Phys. Rev. B 81, 075314 (2010).
  * Carmele  _et al._ (2009) A. Carmele, A. Knorr, and M. Richter, Physical Review B 79, 035316 (2009).
  * Ramsay  _et al._ (2010b) A. J. Ramsay, A. V. Gopal, E. M. Gauger, A. Nazir, B. W. Lovett, A. M. Fox, and M. S. Skolnick, Phys. Rev. Lett. 104, 017402 (2010b).
  * Glässl  _et al._ (2011b) M. Glässl, M. D. Croitoru, A. Vagov, V. M. Axt, and T. Kuhn, Phys. Rev. B 84, 125304 (2011b).
  * Machnikowski and Jacak (2004) P. Machnikowski and L. Jacak, Phys. Rev. B 69, 193302 (2004).
  * Krügel  _et al._ (2006) A. Krügel, V. M. Axt, and T. Kuhn, Phys. Rev. B 73, 035302 (2006).
  * Reiter (2017) D. E. Reiter, Phys. Rev. B 95, 125308 (2017).
  * Naumann  _et al._ (2016) N. L. Naumann, L. Droenner, W. W. Chow, J. Kabuss, and A. Carmele, JOSA B 33, 1492 (2016).
  * Droenner  _et al._ (2017) L. Droenner, N. L. Naumann, J. Kabuss, and A. Carmele, Physical Review A 96, 043805 (2017).
  * Kabuss  _et al._ (2012) J. Kabuss, A. Carmele, T. Brandes, and A. Knorr, Phys. Rev. Lett. 109, 054301 (2012).
  * Kepesidis  _et al._ (2013) K. V. Kepesidis, S. D. Bennett, S. Portolan, M. D. Lukin, and P. Rabl, Phys. Rev. B 88, 064105 (2013).
  * Czerniuk  _et al._ (2017) T. Czerniuk, D. Wigger, A. V. Akimov, C. Schneider, M. Kamp, S. Höfling, D. R. Yakovlev, T. Kuhn, D. E. Reiter, and M. Bayer, Phys. Rev. Lett. 118, 133901 (2017).
  * Richter and Knorr (2010) M. Richter and A. Knorr, Annals of Physics 325, 711 (2010).
  * Stauber  _et al._ (2000) T. Stauber, R. Zimmermann, and H. Castella, Phys. Rev. B 62, 7336 (2000).
  * Heitz  _et al._ (1999) R. Heitz, I. Mukhametzhanov, O. Stier, A. Madhukar, and D. Bimberg, Phys. Rev. Lett. 83, 4654 (1999).
  * Heitz  _et al._ (2001) R. Heitz, H. Born, F. Guffarth, O. Stier, A. Schliwa, A. Hoffmann, and D. Bimberg, Phys. Rev. B 64, 241305 (2001).
  * Droenner  _et al._ (2019) L. Droenner, R. Finsterhölzl, M. Heyl, and A. Carmele, arXiv preprint arXiv:1902.04986 (2019).
  * Leggett  _et al._ (1987) A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
  * Vagov  _et al._ (2006) A. Vagov, M. Croitoru, V. Axt, T. Kuhn, and F. Peeters, physica status solidi (b) 243, 2233 (2006).
  * Vagov  _et al._ (2011) A. Vagov, M. Croitoru, V. Axt, P. Machnikowski, and T. Kuhn, physica status solidi (b) 248, 839 (2011).
  * Dousse  _et al._ (2008) A. Dousse, L. Lanco, J. Suffczyński, E. Semenova, A. Miard, A. Lemaître, I. Sagnes, C. Roblin, J. Bloch, and P. Senellart, Phys. Rev. Lett. 101, 267404 (2008).
  * Gschrey  _et al._ (2013) M. Gschrey, F. Gericke, A. Schüßler, R. Schmidt, J.-H. Schulze, T. Heindel, S. Rodt, A. Strittmatter, and S. Reitzenstein, Applied Physics Letters 102, 251113 (2013).
  * Sotier  _et al._ (2009) F. Sotier, T. Thomay, T. Hanke, J. Korger, S. Mahapatra, A. Frey, K. Brunner, R. Bratschitsch, and A. Leitenstorfer, Nature Physics 5, 352 (2009).
  * Flagg  _et al._ (2009) E. B. Flagg, A. Muller, J. Robertson, S. Founta, D. Deppe, M. Xiao, W. Ma, G. Salamo, and C.-K. Shih, Nature Physics 5, 203 (2009).
  * Laussy  _et al._ (2008) F. P. Laussy, E. Del Valle, and C. Tejedor, Physical review letters 101, 083601 (2008).
  * Kasprzak  _et al._ (2010) J. Kasprzak, S. Reitzenstein, E. A. Muljarov, C. Kistner, C. Schneider, M. Strauss, S. Höfling, A. Forchel, and W. Langbein, Nature materials 9, 304 (2010).
  * Kavokin (2017) A. Kavokin,  _Microcavities_ , 21 (Oxford University Press, 2017).
  * Lodahl  _et al._ (2015) P. Lodahl, S. Mahmoodian, and S. Stobbe, Rev. Mod. Phys. 87, 347 (2015).
  * Schneebeli  _et al._ (2008) L. Schneebeli, M. Kira, and S. W. Koch, Phys. Rev. Lett. 101, 097401 (2008).
  * Hopfmann  _et al._ (2017) C. Hopfmann, A. Carmele, A. Musiał, C. Schneider, M. Kamp, S. Höfling, A. Knorr, and S. Reitzenstein, Physical Review B 95, 035302 (2017).
  * Lounis and Orrit (2005) B. Lounis and M. Orrit, Reports on Progress in Physics 68, 1129 (2005).
  * Ding  _et al._ (2016) X. Ding, Y. He, Z.-C. Duan, N. Gregersen, M.-C. Chen, S. Unsleber, S. Maier, C. Schneider, M. Kamp, S. Höfling, C.-Y. Lu, and J.-W. Pan, Phys. Rev. Lett. 116, 020401 (2016).
  * Somaschi  _et al._ (2016) N. Somaschi, V. Giesz, L. De Santis, J. Loredo, M. P. Almeida, G. Hornecker, S. L. Portalupi, T. Grange, C. Antón, J. Demory, _et al._ , Nature Photonics 10, 340 (2016).
  * Iles-Smith  _et al._ (2017) J. Iles-Smith, D. P. McCutcheon, A. Nazir, and J. Mørk, Nature Photonics 11, 521 (2017).
  * Calic  _et al._ (2011) M. Calic, P. Gallo, M. Felici, K. Atlasov, B. Dwir, A. Rudra, G. Biasiol, L. Sorba, G. Tarel, V. Savona, _et al._ , Physical review letters 106, 227402 (2011).
  * Kaer  _et al._ (2010) P. Kaer, T. R. Nielsen, P. Lodahl, A.-P. Jauho, and J. Mørk, Phys. Rev. Lett. 104, 157401 (2010).
  * Hughes  _et al._ (2011a) S. Hughes, P. Yao, F. Milde, A. Knorr, D. Dalacu, K. Mnaymneh, V. Sazonova, P. J. Poole, G. C. Aers, J. Lapointe, R. Cheriton, and R. L. Williams, Phys. Rev. B 83, 165313 (2011a).
  * Winger  _et al._ (2009) M. Winger, T. Volz, G. Tarel, S. Portolan, A. Badolato, K. J. Hennessy, E. L. Hu, A. Beveratos, J. Finley, V. Savona, _et al._ , Physical review letters 103, 207403 (2009).
  * Hughes  _et al._ (2011b) S. Hughes, P. Yao, F. Milde, A. Knorr, D. Dalacu, K. Mnaymneh, V. Sazonova, P. Poole, G. Aers, J. Lapointe, _et al._ , Physical Review B 83, 165313 (2011b).
  * Carmele  _et al._ (2013b) A. Carmele, J. Kabuss, and W. W. Chow, Physical Review B 87, 041305 (2013b).
  * Thoma  _et al._ (2016) A. Thoma, P. Schnauber, M. Gschrey, M. Seifried, J. Wolters, J.-H. Schulze, A. Strittmatter, S. Rodt, A. Carmele, A. Knorr, T. Heindel, and S. Reitzenstein, Phys. Rev. Lett. 116, 033601 (2016).
  * Rempe  _et al._ (1987) G. Rempe, H. Walther, and N. Klein, Phys. Rev. Lett. 58, 353 (1987).
  * Eberly  _et al._ (1980) J. H. Eberly, N. B. Narozhny, and J. J. Sanchez-Mondragon, Phys. Rev. Lett. 44, 1323 (1980).
  * Hong  _et al._ (1987) C. K. Hong, Z. Y. Ou, and L. Mandel, Phys. Rev. Lett. 59, 2044 (1987).
  * Galperin  _et al._ (2006) Y. M. Galperin, B. L. Altshuler, J. Bergli, and D. V. Shantsev, Phys. Rev. Lett. 96, 097009 (2006).
  * Laikhtman (1985) B. D. Laikhtman, Phys. Rev. B 31, 490 (1985).
  * Eberly  _et al._ (1984) J. H. Eberly, K. Wódkiewicz, and B. W. Shore, Phys. Rev. A 30, 2381 (1984).
  * Liu  _et al._ (2018b) J. Liu, K. Konthasinghe, M. Davanço, J. Lawall, V. Anant, V. Verma, R. Mirin, S. W. Nam, J. D. Song, B. Ma, Z. S. Chen, H. Q. Ni, Z. C. Niu, and K. Srinivasan, Phys. Rev. Applied 9, 064019 (2018b).
  * Glässl  _et al._ (2012) M. Glässl, L. Sörgel, A. Vagov, M. D. Croitoru, T. Kuhn, and V. M. Axt, Phys. Rev. B 86, 035319 (2012).