# Context-Aware Generative Models for Prediction of Aircraft Ground Tracks

Nick Pepper1, George De Ath2, Marc Thomas3, Richard Everson2 and Tim Dodwell2
1 The Alan Turing Institute 2 University of Exeter 3 NATS. Manuscript received
xxxx; revised xxxx. This work was supported by an EPSRC Prosperity Partnership
between NATS, Turing, Exeter and Cambridge. Nick Pepper
The Alan Turing Institute
The British Library
London, UK
<EMAIL_ADDRESS>
&George De Ath
Department of Computer Science
University of Exeter
Exeter, UK
&Marc Thomas
NATS
Whitely, Fareham, UK
&Richard Everson
Department of Computer Science
University of Exeter
Exeter, UK
&Tim Dodwell
Department of Engineering
University of Exeter
Exeter, UK

###### Abstract

Trajectory prediction (TP) plays an important role in supporting the decision-
making of Air Traffic Controllers (ATCOs). Traditional TP methods are
deterministic and physics-based, with parameters that are calibrated using
aircraft surveillance data harvested across the world. These models are,
therefore, agnostic to the intentions of the pilots and ATCOs, which can have
a significant effect on the observed trajectory, particularly in the lateral
plane. This work proposes a generative method for lateral TP, using
probabilistic machine learning to model the effect of the epistemic
uncertainty arising from the unknown effect of pilot behaviour and ATCO
intentions. The models are trained to be specific to a particular sector,
allowing local procedures such as coordinated entry and exit points to be
modelled. A dataset comprising a week’s worth of aircraft surveillance data,
passing through a busy sector of the United Kingdom’s upper airspace, was used
to train and test the models. Specifically, a piecewise linear model was used
as a functional, low-dimensional representation of the ground tracks, with its
control points determined by a generative model conditioned on partial
context. It was found that, of the investigated models, a Bayesian Neural
Network using the Laplace approximation was able to generate the most
plausible trajectories in order to emulate the flow of traffic through the
sector.

##  1 Introduction

Air Traffic Control (ATC) issues instructions to aircraft in order to prevent
collisions by ensuring adequate separation between aircraft, as well as
enabling the expeditious and orderly flow of air traffic [1]. In order to
achieve this, the airspace is divided into multiple interlocking polyhedra of
differing sizes, known as sectors, with the aircraft in each sector usually
controlled by a single Air Traffic Control Officer (ATCO) with the support of
a planning controller.

In order to meet future requirements, extensive ATC modernisation programmes
[2, 3] are ongoing, which seek to develop a number of advanced tools capable,
for example, of providing decision support, improving workload forecasting and
planning, enabling more fuel efficient procedures and systems, optimising
route networks, improving network traffic predictions, and providing more
accurate conflict detection tools [4].

A prerequisite for many of these advanced tools is an accurate Trajectory
Prediction (TP), indicating the expected path an aircraft will take.
Traditional TP models used within ATC are physics-based (e.g. [5]), using
differential equations describing flight mechanics in order to estimate the
aircraft’s future positions. However, these deterministic physics-based models
are insufficient for all the requirements of the modernisation initiatives, as
they do not model the ATCO instructions or specific airline procedures, and,
therefore, have limited utility in longer distance trajectories if these are
not explicitly known. Current research efforts focus on TP methods that can
implicitly model these factors.

###  1.1 Related work

Several previous papers have used data-driven methods to predict complete
trajectories, often from the origin to the destination airports, implicitly
modelling the planned routes, ATCO instructions, and procedures followed.
Examples of this include the SESAR-funded DART project [6, 7], which examined
a range of methods including Hidden Markov Models [8] to model origin to
destination trajectories without explicitly modelling ATCO instructions.
Functional representations of aircraft trajectories, derived from real-world
data, have also been explored for detecting outliers [9] and for forecasting
the future trajectories of aircraft [10, 11].

Increasingly, research has focused on developing machine learning methods for
the task of 4D trajectory prediction in the short- to medium-term, leveraging
the large quantity of radar surveillance data available to ATC [12, 13, 14,
15, 16, 17]. However, this is a challenging task due to the effects of weather
uncertainties [18, 19, 20], changing intentions of ATCOs, and the effect of
local procedures [21]. To date, such models have simplified this complexity by
training machine learning models to emulate aircraft following a single route,
or a limited set of routes. One solution to this issue is to use a hybrid
method to ‘evolve’ the deterministic, physics-based models that are currently
used in TP. For instance, Hratsovec and Solina [22] use a $k$-nearest
neighbour classifier to account for the effect of airline procedures by
conditioning the parameters used in BADA on the filed flight plan. Similarly,
a multiple model approach has been suggested in other works (see, e.g. [23])
for modelling changing pilot intent. The generative model proposed in this
work follows a similar approach by modelling the ground track followed by the
aircraft, which is largely conditioned on intent, separately from the aircraft
speed and climb rate, which are also affected by meteorological conditions and
subject to constraints imposed by flight mechanics.

###  1.2 Contribution

This work focuses on developing a data-driven generative model to replicate
the lateral motion of aircraft through a specific sector, employing machine
learning techniques to leverage the large datasets of aircraft surveillance
data available to ATC, along with the planned route of the aircraft, which are
filed prior to take-off. This enables prediction of the path of aircraft
through this sector, implicitly taking into account the air traffic control
instructions issued and sector- or route-specific procedures followed, which
would otherwise be challenging to model explicitly. The generative nature of
the model allows a range of realistic trajectories to be produced. By coupling
these generated trajectories with existing, physics-based aircraft models such
as BADA for a full 4D trajectory prediction, the model could be useful, for
example, in tools for optimising arrival times, improving network flows or
predicting potential traffic hotspots.

This work evaluates the suitability of Bayesian Neural Networks (BNNs) [24,
25, 26, 27] and Deep Ensembles (DEs) [28] for this task and improves upon the
works discussed in the previous section in several respects:

  1. 1.

Trajectories are generated in the lateral plane that are conditioned on the
filed flight plan and implicitly take into account ATCO intentions

  2. 2.

A probabilistic machine learning model is trained for multiple routes in one
of the busiest sectors of airspace in the UK, including converging and
diverging routes

  3. 3.

As a generative model, the method is assessed by how well it can reconstruct
the flow of traffic through a sector, rather than on a trajectory-by-
trajectory basis

  4. 4.

The method is modular and can be easily coupled with existing physics-based
models, such as BADA, to give a full 4D prediction

In the next section the piecewise linear representation used to model aircraft
ground tracks is introduced and the BNN and DE architectures are outlined.
Section 3 describes the dataset of real-world surveillance data that was used
to train and test the model. In Section 4 a pair of metrics for assessing the
quality of the probabilistic trajectory predictions is outlined. The
performance of various probabilistic machine learning methods on the test
dataset are then compared, using these metrics.

##  2 Methods

The piecewise linear representation used to model aircraft ground tracks is
introduced in Section 2.1. Section 2.2 outlines the Bayesian Neural Network
and Deep Ensemble methods that were used as generative models in the piecewise
linear representation.

###  2.1 Piecewise linear representation of aircraft ground tracks

Let $\bm{x}(t)\subseteq\mathbb{X}\in\Re^{2}$ denote an aircraft’s position, at
time $t$, in Geographical Coordinate Space (GCS) that passes through a sector
of airspace, denoted $\Omega\subseteq\mathbb{X}$. The ground track of a
trajectory in this airspace can be discretised as the sequence of radar
observations $X=[\bm{x}^{(0)},\,\bm{x}^{(1)},\dots,\bm{x}^{(n)}]$, observed at
times $\tau_{0}=0,\,\tau_{1},\dots,\tau_{n}$. The distance between two points
in GCS is computed using the Haversine distance, $\mathcal{D}(\bm{a},\bm{b})$
with $\bm{a},\bm{b}\in\mathbb{X}$.

Exploiting the spatio-temporal correlation of the points in these
trajectories, we define a piecewise linear basis to provide a compact
representation of the trajectory. In this representation, ground tracks are
piecewise linear paths that connect a set of control points
$P\in[\bm{p}_{0},\bm{p}_{1},\dots\bm{p}_{d}]$, where
$\bm{p}_{i}\subseteq{\mathbb{X}}$ refers to the location of the $i$-th turn.
The maximum number of turns taken, $d$, is referred to as the degree of the
piecewise representation. Similarly to the control points of Bézier curves
[29], $\bm{p}_{0}$ is defined to be the first point of the trajectory, i.e.,

$\displaystyle\bm{p}_{0}:=\bm{x}^{(0)}.$ (1)

To account for the unequal lengths of the trajectories in the training data,
each trajectory is interpolated onto a normalised timescale $t\in[0,1]$, with:

$\displaystyle\bm{x}(t_{0})=\bm{p}_{0}\;\text{and}\;\bm{x}(t_{d})=\bm{p}_{d},$
(2)

where $t_{0}=0$ and $t_{d}=1$. Arrival times $t_{i},\,i=1,\dots d-1$ at the
control points collected in $P$, using the normalised timescale, are
proportional to the length of the trajectory travelled:

$\displaystyle
t_{i}=t_{i-1}+\frac{\mathcal{D}(\bm{p}_{i},\bm{p}_{i-1})}{\sum_{i=1}^{d}\mathcal{D}(\bm{p}_{i},\bm{p}_{i-1})},\;i=1:d-1.$
(3)

Using the arrival times and control points in $P$, the GCS coordinates of a
point in the piecewise linear representation can be determined through:

$\displaystyle\bm{x}(t)=\bm{x}(t_{i})+\frac{t-t_{i}}{t_{i+1}-t_{i}}({\bm{p}}_{i+1}-{\bm{p}}_{i}),$
(4)

where ${t}_{i}\leq t\leq{t}_{i+1}$, with $t_{i}$ referring to the time at
which the last control point in $P$ was reached and $t_{i+1}$ to the arrival
time at the next control point. For an individual trajectory, $X$, in the
training data, an optimisation process was followed to determine the optimal
locations of the points in $P$ that minimised a cost function
$\mathcal{J}(\cdot)$:

$\displaystyle\hat{P}=\underset{P}{\textrm{arg min}}\;\mathcal{J}(P|X).$ (5)

$\mathcal{J}(\cdot)$ quantifies how well the piecewise linear representation
fits the radar data in $X$, and is defined as:

$\displaystyle\mathcal{J}(P)=\sum_{j=0}^{n}\Big{|}\Big{|}\bm{x}^{(j)}-\hat{\bm{x}}^{(j)}\Big{|}\Big{|}^{2}+\lambda\mathcal{I}(\hat{X}),$
(6)

where $\hat{\bm{x}}^{(j)}$ refers to a point at time $t=\tau^{(j)}/\tau^{(n)}$
along the trajectory generated by (3) and (4), conditioned on $P$, and
collected as $\hat{X}$. The cost function consists of two components, a data
fit term in the form of the $L2$ distance between $X$ and $\hat{X}$, and a
penalty term to dissuade trajectories from making many large turns. The
indicator function, $\mathcal{I}$ is defined as:

$\displaystyle\mathcal{I}(\hat{X})=\begin{cases}0\;\;\text{if}\;\phi_{\scaleobj{0.75}{\Sigma}}(\hat{X})\leq\phi_{u}\\\
1\;\;\text{otherwise}\end{cases}\\!\\!\\!\\!,$ (7)

where $\phi_{\scaleobj{0.75}{\Sigma}}(\hat{X})$ represents the total angle
turned through by the trajectory and $\phi_{u}$ is a user-selected maximum
bound on this quantity. The regularisation constant, $\lambda$, is set to a
value large enough that if the total angle turned by the trajectory exceeds
$\phi_{u}$ a severe penalty is given. $\lambda$ was adjusted to be an order of
magnitude greater than the typical $L2$ distance between $X$ and $\hat{X}$ in
the training data. Similarly, $\phi_{u}$ was set to the largest
$\phi_{\scaleobj{0.75}{\Sigma}}$ observed in the training data.

###  2.2 Generative model for aircraft ground tracks

This section presents a sector-specific generative model for aircraft ground
tracks based on an aircraft’s filed flight plan and its state upon entering.
This information is collected in the vector $\bm{\xi}$. The model generates
ground tracks by performing the probabilistic mapping
$\bm{\xi}\rightarrow\\{\bm{p_{i}}\\}_{i=1}^{d}$, where
$P=\\{\bm{p_{i}}\\}_{i=0}^{d}$ contains a sequence of control points in the
piecewise linear representation discussed above. Here, $\bm{p}_{0}$ represents
the GCS coordinates at which the aircraft enters $\Omega$ and anchors the
generated ground track to the aircraft’s initial location. The conditional
distribution $P|\bm{\xi}$ can be used to generate the ground track using (3)
and (4). The quantities collected in $\bm{\xi}$ are:

  * •

The aircraft position $\bm{p}_{0}$ in GCS coordinates on entering ${\Omega}$

  * •

The aircraft flight level (altitude) when entering $\Omega$

  * •

The filed flight plan (encoded as a sequence of GCS coordinates)

$\bm{\xi}$ necessarily only contains partial information about the context of
the trajectory. The strategy of the ATCO for handling the traffic in their
sector will also affect the trajectory of the aircraft. This will, in part, be
determined by factors such as the time of day, the number of aircraft in the
sector, and the state of neighbouring sectors. It is, therefore, difficult to
encode the strategy of the ATCO directly in the model. For this reason,
probabilistic machine learning is used to inject uncertainty into the model
predictions that accounts for this epistemic uncertainty. In this paper the
suitability of Deep Ensembles [28] and Bayesian Neural Networks [24, 25, 26,
27] for this task are investigated. Both methods are popular ways of adding
uncertainty to neural network predictions in the literature [30]. We outline
the main features of the methods below. For convenience, we denote the vector
that collects the set $\\{\bm{p_{i}}\\}_{i=1}^{d}$ as $\bm{y}\in\Re^{2d}$.
Evaluating (6) for every trajectory in a dataset of $n_{f}$ individual flights
through $\Omega$ yields the training set
$D=\\{\bm{\xi}^{(k)},\bm{y}^{(k)}\\},\,k=1:n_{f}$.

Deep Ensembles (DEs)
Lakshminarayanan et al. [28] proposed to use an ensemble of $n_{e}$
independently-trained neural networks as a way to estimate predictive
uncertainty using neural networks. Each of the $n_{e}$ networks consist of the
same architecture, but are initialised with different weights
$\bm{\theta}\in\Re^{n_{{\theta}}}$. Once trained, the mean and standard
deviation of the networks’ outputs are used as the parameters of i.i.d.
Gaussian distributions that represent the ensemble’s average prediction and
associated uncertainty. In this work, we model each component of $\bm{y}$ as
the output of the deep ensemble. Each network is trained independently to
minimise the loss function:

$\displaystyle\hat{\bm{\theta}}=\underset{\bm{\theta}}{\text{arg
min}}\;\mathcal{L}(D|\bm{\theta}),$ (8)

where $\mathcal{L}(\cdot)$ is defined as:

$\displaystyle{\mathcal{L}(D|\bm{\theta})=r(\bm{\theta})+\sum_{k=1}^{n_{f}}l\big{(}f(\bm{\xi}^{(k)}|\bm{\theta}),\bm{y}^{(k)}\big{)}},$
(9)

and where $r(\cdot)$ is a regularisation term, $f(\cdot)$ represents the
output of the network, and $l(\cdot)$ is the negative-log likelihood. Having
trained the $n_{e}$ networks in the ensemble, using the Adam optimiser [31],
the predictions of the ensemble for an unseen context, $\bm{\xi}$, are
represented as the set of i.i.d. Gaussian probability distributions
$\bm{y}_{m}=N(\mu^{*}_{m},\sigma^{*}_{m}),\,m=1:2d$, with:

$\displaystyle\mu^{*}_{m}(\bm{\xi})$
$\displaystyle=\frac{1}{n_{e}}\sum_{o=1}^{n_{e}}\mu_{o,m}(\bm{\xi}),$ (10)
$\displaystyle\sigma^{*2}_{m}(\bm{\xi})$
$\displaystyle={\sigma^{2}_{\mu_{m}}(\bm{\xi})+\sigma^{2}_{\sigma_{m}}(\bm{\xi})},$

where the standard deviation of the ensemble contains contributions from the
variance of the individual networks, $\sigma^{*}_{\sigma}$, and the spread of
their means, $\sigma^{*}_{\mu}$, with:

$\displaystyle\sigma^{2}_{\mu_{m}}(\bm{\xi})$
$\displaystyle=\frac{1}{n_{e}}\sum_{o=1}^{n_{e}}\mu^{2}_{o,m}(\bm{\xi})-\mu_{m}^{*2}(\bm{\xi}),$
(11) $\displaystyle\sigma^{2}_{\sigma_{m}}(\bm{\xi})$
$\displaystyle=\frac{1}{n_{e}}\sum_{o=1}^{n_{e}}\sigma_{o,m}^{2}(\bm{\xi}).$

The ensemble can be used to generate trajectories by sampling the coordinates
of the control points from the i.i.d. Gaussian distributions, and using (3)
and (4) to determine the ground track that connects these control points.

Bayesian Neural Networks (BNNs)
Deep ensembles incorporate uncertainty into neural network predictions by
expressing the outputs as Gaussian distributions, rather than as point
estimates. A drawback of this approach is that each element in posterior
samples of $\bm{y}$ are sampled independently. An alternative approach is to
add uncertainty by making all [27] or a subset [32] of the network parameters,
$\bm{\theta}$, probabilistic. Bayesian methods are then used to infer these
probability distributions, based on the evidence of the training data. The
advantage of this formulation is that it provides a richer probabilistic model
for $\bm{y}$, because the outputs are no longer constrained to be i.i.d.
Gaussian distributions. On the other hand, BNN training is more complex
compared to deep ensembles because standard backpropagation methods cannot be
applied [33]. Methods such as variational inference are commonly criticised
for being challenging to implement and expensive to train [34]. For this
reason, in this work, the posteriors of the uncertain network parameters were
estimated using the Laplace approximation [32]. This is a two-step process:
first a _maximum a posteriori_ (MAP) estimate is generated for $\bm{\theta}$
by solving (8) and (9), with the $L2$ distance used for $l(\cdot)$. Having
obtained the MAP estimate, denoted $\bm{\theta}_{MAP}$, the posterior of the
uncertain parameters is locally approximated with a Gaussian distribution.
This distribution is centred at $\bm{\theta}_{MAP}$, with its covariance
matrix corresponding to the local curvature:

$\displaystyle
p(\bm{\theta}|D)=N(\bm{\theta}|\bm{\theta}_{MAP},\Sigma)\;\text{with}\;\Sigma:=(\nabla_{\bm{\theta}}^{2}\mathcal{L}(D|\bm{\theta})_{\bm{\theta}_{MAP}})^{-1},$
(12)

and where $p(\cdot)$ represents the posterior probability of the network
parameters. The Laplace package [32] was used to implement this post-hoc step,
using a zero-mean Gaussian prior for the uncertain network parameters. The
generalized Gauss-Newton matrix (GGN) was computed to store the derivatives of
the network parameters. Two strategies were deployed to reduce the
computational cost of this step: firstly, by only considering a subset of the
network parameters to be uncertain and secondly, by employing factorization
assumptions for the GGN matrix. Following an initial exploratory phase, the
best results were found using a diagonal factorization with a last-layer
Laplace approximation.

Linear model
A probabilistic linear model is used as a baseline with which to compare the
machine learning methods to. For lateral TP, a first assumption is that
aircraft will fly directly to the exit of the sector. In the linear model used
here, this exit point is the last waypoint in the flight plan. To make the
model probabilistic, jitter was added to the GCS coordinates of the exit
point, using a two-dimensional, zero-mean, Gaussian distribution with
covariance $\Sigma=\sigma^{2}I$, where $\sigma=0.05$°.

##  3 Data Preparation

Probabilistic machine learning models were trained to emulate the flow of air
traffic through a specific sector located in the upper airspace of the UK,
above several of the world’s busiest airports. $\Omega$ was defined to include
the sector and a 50 nautical mile buffer around it’s edge. The sector boundary
in the lateral plane is illustrated in Fig. 1. The dataset contained one
week’s worth of radar surveillance data, filtered to the 1902 aircraft that
were under the control of ATCOs assigned to this sector during this period.
For ease of comparison, another filtering step was performed: selecting those
flights that followed the eight most popular filed flight plans in the
dataset. Of the 1152 remaining trajectories, an 80:20 training:test split was
performed, resulting in a training dataset of 922 flights and a test dataset
of 230 flights. The top-left panel of Fig. 1 displays the GCS coordinates of
the waypoints in the filed routes, together with lines connecting the
waypoints of each unique flight plan.

The ground tracks of these trajectories are illustrated in the lower panels of
Fig. 1, with lines coloured by the filed flight plan. As can be seen from Fig.
1, aircraft are not constrained to fly from one fix to the next and there can
be considerable variation between the route joining together the waypoints and
the actual trajectory followed. Many of these routes are overlapping, with
some routes sharing waypoints. Flows of traffic in the sector are dominated by
aircraft either departing from or arriving at major airports in the south of
the UK. The top-right panel of Fig. 1 indicates the flight level (altitude) of
the aircraft on the normalised timescale. Routes 1 and 8 correspond to
descending aircraft entering the sector from the east and south, respectively,
and heading north. Routes 2-7 correspond to climbing aircraft that enter the
sector from the north, having departed from these airports, and are heading
south.

The trajectories in the training dataset were projected onto the piecewise
linear basis discussed in Section 2, with the control points determined by
solving (6) for each trajectory. Based on an inspection of the typical number
of turns taken by aircraft in this dataset, the degree of the representation
was set to $d=3$.

Figure 1: Boundary of the studied sector, locations of the waypoints, and the
eight unique routes through $\Omega$ in the dataset (top-left). Flight levels
of the trajectories in the dataset, interpolated onto $t\in[0,1]$ (top-right).
Trajectories in the dataset, coloured by their flight plan through $\Omega$
(bottom).

##  4 Application to a sector of UK airspace

The probabilistic machine learning models described in Section 2 were trained
using the training set outlined in the previous section, and, for each of the
test trajectories, $\bm{\xi}$, 10 posterior samples were drawn. In order to
assess how well the methods emulated the flow of traffic, it was necessary to
quantify the difference between the generated trajectories and the
trajectories in the test dataset. While there are numerous methods in the
literature for computing the distance between trajectories (see, e.g. [35]),
there is no fixed convention for measuring trajectory similarity in ATC.
Common metrics for measuring trajectory similarity such as the Euclidean
distance, Hausdorff distance [36], and Dynamic Time Warping [37] implicitly
account for differences in heading. However, in ATC the relative heading
between trajectories is also an important consideration. As an example,
trajectories in routes 3 and 4 overlap for much of the sector, but, from an
ATC perspective, would be considered to be distinct because they exit the
sector across different boundaries. Summarising the similarity of trajectories
with a single statistic would require finding a weighting between the relative
headings of trajectories and their spatial proximity for ATC. For this reason,
we simplify the measurement of trajectory similarity by assessing how well the
probabilistic methods emulate the flow of air traffic through a sector using a
geometric method. The advantage of this approach is that the spatial proximity
and relative headings between two sets of trajectories can be quantified
separately.

Fig. 2 is a schematic that illustrates the method using the training data. For
a given route, in this case route 2, a plane is defined that is perpendicular
to the filed flight plan. The origin of this plane is the last waypoint that
is within the sector boundary. This plane is represented by a dashed blue line
in Fig. 2. The Haversine distance, in nautical miles, between the origin and
the point at which a trajectory intersects this plane is calculated, along
with the heading of the trajectory at this point. These quantities are denoted
$\mathcal{D}_{H}$ and $\phi$ respectively. Histograms of $\mathcal{D}_{H}$ and
$\text{sin}(\phi)$ can then be computed for each of the eight studied routes.
The sine of the heading is used in order to make the headings continuous. The
schematic in the right panel of Fig. 2 illustrates the histogram extracted for
$\mathcal{D}_{H}$ from the training data for route 2. For a given route, the
similarity of two sets of trajectories is assessed by comparing the
statistical distance between the histograms of these two quantities.

Figure 2: Schematic illustrating the geometric method used to assess how well
the proposed generative emulated the distribution of ground tracks in the
dataset.

Fig. 3 illustrates the trajectories generated by the BNN (top) and the Deep
Ensemble (bottom). The plots also show the waypoints of each route as coloured
dots, joined by solid lines. The dashed coloured lines indicate the planes
that were used to compute $\mathcal{D}_{H}$ and $\phi$. The posterior samples
from the Deep Ensemble are more widely distributed across the sector, however,
the consequence of this is that many samples make turns that are not observed
in the surveillance data. Conversely, the posterior samples from the BNN are
more tightly clustered, but, from inspection, appear both more similar to the
trajectories in Fig. 1, and more physically plausible. Fig. 4 displays kernel
density estimates for the Probability Density Functions (PDFs) of
$\mathcal{D}_{H}$ for the test dataset, the BNN, Deep Ensemble, and the
probabilistic linear model. Similarly, Fig. 5 displays kernel density
estimates of the PDFs for $\text{sin}(\phi)$ for the test data and the various
methods. Tables 1 and 2 tabulate the statistical distance between the
probability distributions of the data and the three methods, as quantified by
the Kolmogorov-Smirnov [38] (KS) distance and percentage difference of the
mean respectively. The statistical distance for the best performing model is
in bold for both quantities and each of the eight routes.

As might be expected, the linear method performed best for route 8, where the
majority of trajectories are straight lines across the sector, but it does not
perform well for the other routes which have more complex trajectories. In
general, the distribution of $\mathcal{D}_{H}$ for the Deep Ensemble is closer
to that of the test data than the BNN. The distributions from the BNN tend to
accurately estimate the mean but, except for the first and eighth routes, it
underestimates the variance. However, the headings at the plane intersection
of the trajectories generated by the BNN are much closer to the test dataset
than the other two generative models. This is likely due to the fact that
samples generated from the BNN have correlations between elements in $\bm{y}$.
In contrast, the samples for the $2d$ elements of $\bm{y}$ from the DE are
independent, and, therefore, often introduce spurious turns into the generated
trajectories; c.f. route 8 in Fig. 3.

Figure 3: Generated trajectories, sampled from the posteriors of the Bayesian Neural Network (top) and Deep Ensemble (bottom). Solid black lines indicate the waypoints of the route, with dotted lines indicating the plane on which $\mathcal{D}_{H}$ and $\phi$ are computed for each route. Figure 4: Kernel density plot of $\mathcal{D}_{H}$, coloured by unique route through $\Omega$, for the test dataset and the trajectories generated by the probabilistic models. Table 1: Statistical distance between distributions for $\mathcal{D}_{H}$ and $\text{sin}(\phi)$, as quantified by the KS distance, for each of the eight unique routes.  Route | KS distance ($\mathcal{D}_{H}$) $\downarrow$ | KS distance ($\text{sin}(\phi)$) $\downarrow$
---|---|---
| Lin. | DE | BNN | Lin. | DE | BNN
1 | 1.000 | 0.483 | 0.578 | 1.000 | 0.378 | 0.178
2 | 0.905 | 0.457 | 0.805 | 0.848 | 0.267 | 0.352
3 | 0.648 | 0.188 | 0.472 | 0.772 | 0.344 | 0.356
4 | 0.912 | 0.235 | 0.238 | 0.527 | 0.262 | 0.142
5 | 0.553 | 0.138 | 0.233 | 0.420 | 0.171 | 0.111
6 | 1.000 | 0.200 | 0.464 | 0.966 | 0.205 | 0.150
7 | 0.482 | 0.206 | 0.406 | 0.341 | 0.247 | 0.176
8 | 0.908 | 0.479 | 0.692 | 0.562 | 0.404 | 0.179
Table 2: Statistical distance between distributions for $\mathcal{D}_{H}$ and $\text{sin}(\phi)$, as quantified by the percentage error in the mean, for each of the eight unique routes. Route |  $\Delta\mathbb{E}(\mathcal{D}_{H}$) (%) $\downarrow$ |  $\Delta\mathbb{E}(\text{sin}(\phi))$ (%) $\downarrow$
---|---|---
| Lin. | DE | BNN | Lin. | DE | BNN
1 | 506.45 | 122.89 | 129.61 | 55.43 | 26.16 | 2.7
2 | 67.48 | 30.08 | 35.47 | 10.47 | 18.74 | 0.58
3 | 74.58 | 5.04 | 27.19 | 3.44 | 0.42 | 2.38
4 | 238.51 | 19.65 | 5.11 | 1121.90 | 134.94 | 227.64
5 | 58.20 | 1.06 | 0.37 | 1.03 | 0.02 | 3.64
6 | 289.74 | 27.19 | 8.45 | 13.56 | 9.33 | 5.43
7 | 68.97 | 1.17 | 14.23 | 0.44 | 0.47 | 0.64
8 | 1352.92 | 132.48 | 739.97 | 11.17 | 19.10 | 2.80
Figure 5: Kernel density plot of $\text{sin}(\phi)$, coloured by unique route
through $\Omega$, for the test dataset and the trajectories generated by the
probabilistic models.

##  5 Conclusions

A probabilistic machine learning method has been presented for generating
plausible ground tracks for an aircraft entering a specific sector of
airspace. The model is conditioned on contextual information gathered from
real-world aircraft surveillance data. A piecewise linear representation was
used to model the ground tracks. It was found that, of the probabilistic
models evaluated, the best performance was achieved using a Bayesian Neural
Network, with posterior probability distributions over the weights in the
last-layer provided by the Laplace Approximation. This is likely due to the
natural way in which correlations between the locations of control points can
be encoded within the model. It was found that this probabilistic model could
be used to plausibly emulate the flow of air traffic through a busy sector of
UK airspace with real-world data.

Future work includes coupling the presented probabilistic models for aircraft
ground tracks with a deterministic, physics-based TP model such as the Base of
Aircraft Data (BADA) model [5]. Combined, these models would provide a
probabilistic method for four-dimensional TP (GCS coordinates, time, and
altitude) that would be capable of propagating the effect of epistemic
uncertainty arising from unknown pilot and ATCO intentions.

## References

  * [1] ICAO, “Doc 4444 pans-atm,” _International Civil Aviation Organization_ , 2001\.
  * [2] Federal Aviation Administration. (2023) Next generation air transportation system (nextgen). [Online]. Available: https://www.faa.gov/nextgen
  * [3] SESAR Joint Undertaking. (2021) Delivering the digital european sky. [Online]. Available: https://www.sesarju.eu/
  * [4] ——. (2021) SESAR solutions catalogue (fourth edition). [Online]. Available: https://www.sesarju.eu/activities-solutions
  * [5] A. Nuic, D. Poles, and V. Mouillet, “BADA: An advanced aircraft performance model for present and future ATM systems,” _International Journal of Adaptive Control and Signal Processing_ , vol. 24, no. 10, pp. 850–866, 2010.
  * [6] E. C. Fernández, J. M. Cordero, G. Vouros, N. Pelekis, T. Kravaris, H. Georgiou, G. Fuchs, N. Andrienko, G. Andrienko, E. Casado _et al._ , “DART: A machine-learning approach to trajectory prediction and demand-capacity balancing,” _SESAR Innovation Days, Belgrade_ , pp. 28–30, 2017.
  * [7] SESAR Joint Undertaking. (2018) DART final project results report. [Online]. Available: https://www.sesarju.eu/node/3179
  * [8] L. E. Baum and T. Petrie, “Statistical inference for probabilistic functions of finite state markov chains,” _The Annals of Mathematical Statistics_ , vol. 37, no. 6, pp. 1554–1563, 1966.
  * [9] F. Nicol, “Functional principal component analysis of aircraft trajectories,” in _Proceedings of the 2nd International Conference on Interdisciplinary Science for Innovative Air Traffic Management (ISIATM)_ , 2013.
  * [10] N. Pepper, M. Thomas, G. De Ath, E. Oliver, R. Cannon, R. Everson, and T. Dodwell, “A probabilistic model for aircraft in climb using monotonic functional gaussian process emulators,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 479, no. 2271, pp. 1–19, 2023.
  * [11] K. Tastambekov, S. Puechmorel, D. Delahaye, and C. Rabut, “Aircraft trajectory forecasting using local functional regression in sobolev space,” _Transportation Research Part C: Emerging Technologies_ , vol. 39, pp. 1–22, 2014.
  * [12] S. Ayhan and H. Samet, “Aircraft trajectory prediction made easy with predictive analytics,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , ser. KDD ’16. New York, NY, USA: Association for Computing Machinery, 2016, pp. 21––30.
  * [13] Y. Chen, J. Sun, Y. Lin, G. Gui, and H. Sari, “Hybrid n-inception-lstm-based aircraft coordinate prediction method for secure air traffic,” _IEEE Trans. Intell. Transp. Syst._ , vol. 23, no. 3, pp. 2773–2783, 2022.
  * [14] Z. Shi, M. Xu, and Q. Pan, “4-d flight trajectory prediction with constrained lstm network,” _IEEE Trans. Intell. Transp. Syst._ , vol. 22, no. 11, pp. 7242–7255, 2021.
  * [15] L. Ma and S. Tian, “A hybrid cnn-lstm model for aircraft 4d trajectory prediction,” _IEEE Access_ , vol. 8, pp. 134 668–134 680, 2020.
  * [16] A. De Leege, M. van Paassen, and M. Mulder, “A machine learning approach to trajectory prediction,” in _AIAA Guidance, Navigation, and Control (GNC) Conference_ , 2013, p. 4782.
  * [17] X. Zhang and S. Mahadevan, “Bayesian neural networks for flight trajectory prediction and safety assessment,” _Decision Support Systems_ , vol. 131, p. 113246, 2020.
  * [18] Y. Pang, X. Zhao, H. Yan, and Y. Liu, “Data-driven trajectory prediction with weather uncertainties: A Bayesian deep learning approach,” _Transportation Research Part C: Emerging Technologies_ , vol. 130, p. 103326, 2021.
  * [19] Y. Pang and Y. Liu, “Probabilistic aircraft trajectory prediction considering weather uncertainties using dropout as Bayesian approximate variational inference,” in _AIAA Scitech 2020 Forum_ , 2020, p. 1413.
  * [20] X. Tang, P. Chen, and Y. Zhang, “4d trajectory estimation based on nominal flight profile extraction and airway meteorological forecast revision,” _Aerospace Science and Technology_ , vol. 45, pp. 387–397, 2015.
  * [21] J. Zhang, J. Liu, R. Hu, and H. Zhu, “Online four dimensional trajectory prediction method based on aircraft intent updating,” _Aerospace Science and Technology_ , vol. 77, pp. 774–787, 2018.
  * [22] M. Hrastovec and F. Solina, “Prediction of aircraft performances based on data collected by air traffic control centers,” _Transportation Research Part C: Emerging Technologies_ , vol. 73, pp. 167–182, 2016.
  * [23] J. L. Yepes, I. Hwang, and M. Rotea, “New algorithms for aircraft intent inference and trajectory prediction,” _Journal of guidance, control, and dynamics_ , vol. 30, no. 2, pp. 370–382, 2007.
  * [24] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, “Weight uncertainty in neural networks,” in _Proceedings of the 32nd International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, vol. 37. Lille, France: PMLR, 2015, pp. 1613–1622.
  * [25] Y. Gal, “Uncertainty in deep learning,” Ph.D. dissertation, University of Cambridge, 2016.
  * [26] D. J. C. Mackay, “Probable networks and plausible predictions — a review of practical Bayesian methods for supervised neural networks,” _Network: Computation in Neural Systems_ , vol. 6, no. 3, pp. 469–505, 1995.
  * [27] R. M. Neal, _Bayesian Learning for Neural Networks_. Berlin, Heidelberg: Springer-Verlag, 1996.
  * [28] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles,” in _Advances in Neural Information Processing Systems_ , vol. 30. Curran Associates, Inc., 2017, pp. 6402–6413.
  * [29] G. Farin, _Curves and Surfaces for Computer-Aided Geometric Design_ , 3rd ed. Boston: Academic Press, 1993.
  * [30] A. G. Wilson and P. Izmailov, “Bayesian deep learning and a probabilistic perspective of generalization,” in _Advances in Neural Information Processing Systems_ , vol. 33. Curran Associates, Inc., 2020, pp. 4697–4708.
  * [31] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014.
  * [32] E. Daxberger, A. Kristiadi, A. Immer, R. Eschenhagen, M. Bauer, and P. Hennig, “Laplace redux - effortless Bayesian deep learning,” in _Advances in Neural Information Processing Systems_ , vol. 34. Curran Associates, Inc., 2021, pp. 20 089–20 103.
  * [33] L. V. Jospin, H. Laga, F. Boussaid, W. Buntine, and M. Bennamoun, “Hands-on Bayesian neural networks—a tutorial for deep learning users,” _IEEE Comput. Intell. Mag._ , vol. 17, no. 2, pp. 29–48, 2022.
  * [34] K. Osawa, S. Swaroop, M. E. E. Khan, A. Jain, R. Eschenhagen, R. E. Turner, and R. Yokota, “Practical deep learning with Bayesian principles,” in _Advances in neural information processing systems_ , vol. 32. Curran Associates, Inc., 2019, pp. 4287–4299.
  * [35] Y. Tao, A. Both, R. I. Silveira, K. Buchin, S. Sijben, R. S. Purves, P. Laube, D. Peng, K. Toohey, and M. Duckham, “A comparative analysis of trajectory similarity measures,” _GIScience & Remote Sensing_, vol. 58, no. 5, pp. 643–669, 2021.
  * [36] A. A. Taha and A. Hanbury, “An efficient algorithm for calculating the exact Hausdorff distance,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 37, no. 11, pp. 2153–2163, 2015.
  * [37] D. J. Berndt and J. Clifford, “Using dynamic time warping to find patterns in time series,” in _Proceedings of the 3rd International Conference on Knowledge Discovery and Data Mining_. AAAI Press, 1994, pp. 359–370.
  * [38] F. J. Massey, “The Kolmogorov-Smirnov test for goodness of fit,” _Journal of the American Statistical Association_ , vol. 46, no. 253, pp. 68–78, 1951.


# Non-invasive in silico determination of ventricular wall pre-straining and
characteristic cavity pressures

Sebastian Skatulla  Computational Continuum Mechanics Research Group,
Department of Civil Engineering, University of Cape Town, South Africa  Carlo
Sansour  Department of Mathematics, Bethlehem University, Bethlehem, Palestine
Mary Familusi  Computational Continuum Mechanics Research Group, Department of
Civil Engineering, University of Cape Town, South Africa  South African DST-
NRF Centre of Excellence in Epidemiological Modelling and Analysis,
Stellenbosch University, South Africa  Jagir Hussan  Auckland Bioengineering
Institute, University of Auckland, New Zealand  Ntobeko Ntusi  Division of
Cardiology, Department of Medicine, University of Cape Town, Cape Town, South
Africa  Cape Heart Institute, Faculty of Health Sciences, University of Cape
Town, South Africa  Cape Universities Body Imaging Centre, Faculty of Health
Sciences, University of Cape Town, South Africa  South African Medical
Research Council Extramural Unit on Intersection of Noncommunicable Diseases
and Infectious Diseases

###### Abstract

The clinical application of patient-specific modelling of the heart can
provide valuable insights in supplementing and advancing methods of diagnosis
as well as helping to devise the best possible therapeutic approach for each
individual pathological heart condition. The potential of computational
cardiac mechanics, however, has not yet been fully leveraged due to the
heart’s complex physiology and limitations in the non-invasive in vivo
characterisation of heart properties necessary required for accurate patient-
specific modelling such as the heart anatomy in an unloaded state, ventricular
pressure, the elastic constitutive parameters and the myocardial muscle fibre
orientation distribution. From a solid mechanics point of view without prior
knowledge of the unloaded heart configuration and the cavity pressure-volume
evolution, in particular, the constitutive parameters cannot be accurately
estimated to describe the highly nonlinear elastic material behaviour of
myocardial tissue. Here, knowledge of the volume-normalized end-diastolic
pressure relation for larger mammals is exploited in combination with a novel
iterative inverse parameter optimisation framework to determine end-systolic
and end diastolic pressures, ventricular wall pre-straining and pre-stressing
due the residual end-systolic cavity pressure as well as myocardial tissue
stiffness parameters for biventricular heart models.

##  1 Introduction

Realistic modelling of the whole beating heart is the ”holy grail” of
computational cardiac mechanics which has been achieved only recently by the
international “Living Heart Project” research initiative Baillargeon et al.
(2014). The objective of patient-specific modelling of the heart, however, is
still out of reach, because further intense research efforts are necessary to
determine the biomechanical properties of healthy and diseased heart muscle
tissue in-vivo, meaning non-invasive. The effective implementation and
evaluation of the clinical use of patient-specific modelling supplementing
personalized therapy is mostly unknown and under debate Gray and Pathmanathan
(2018). In patient-specific heart modelling, it is crucial to be able to
calibrate the various model parameters to clinically observed data, because
models can only be relied on and be used for prediction if an observed
physiological state can be represented.

Patient-specific biomechanical modelling of the heart has to account for the
presence of a physiological pressure load as a result of the pre-stressed
state of the imaged tissue. Two main sources of pre-stressing need to be
considered: Firstly, cardiac tissue exhibits growth- and remodelling-related
residual stresses even if completely unloaded ex vivo Omens et al. (2003); Shi
et al. (2019); Zhuan and Luo (2022) and secondly, anatomical heart models are
usually based on in vivo imaging of the beating heart when exhibiting the
lowest cavity pressures at end systole, see e.g. Nikou et al. (2016). The
latter implies that the heart tissue is still experiencing stress, strain and
deformation due to the remaining pressure. Failure to account for this pre-
stressed state in solid tissue mechanics models results in inaccurate metrics,
which are then used for health evaluation, risk assessment, or surgical
planning.

In solid mechanics, the load application and corresponding deformation, strain
and stress response is analyzed as linked to the unloaded material
configuration of the problem at hand. This requires knowledge of the unloaded
configuration which is in cardiac mechanics typically unknown. However, with
the knowledge of the residual cavity pressure, the unknown unloaded
configuration of the heart can be found by solving an inverse problem based on
its known pre-stressed configuration. This inverse problem can usually only be
solved computationally, in particular for biological tissue mechanics
applications due to the highly nonlinear nature. Rajagopal et al. (2007)
proposed a incremental and iterative method applying a finite difference
approximation of the finite element (FE) residual vector formulated in the
Lagrangian form to predicted the undeformed shape of human breasts. Bols et
al. (2013) proposed backward displacement method making use of a fixed point
algorithm to iteratively retrieve the in vivo stress in a mouse-specific
abdominal aorta due to the blood pressure at the moment of imaging and to
reinstate the zero-pressure blood vessel geometry. Using in-vivo pressure data
and magnetic resonance imaging, Pourmodheji et al. (2021) applied a patient-
specific inverse FE modelling framework to iterative determine besides of the
elastic material parameters the pre-stretch of elastin and collagen fibres in
the proximal pulmonary arteries in order to assess mechanical and structural
alterations of these micro-constituents. In contrast, a direct approach
requires the governing equations to be consistently related to the known
deformed configuration. Accordingly, Govindjee and Mihalic (1996) and Gee et
al. (2010) derived the corresponding weak form as integrals over the current
body. Govindjee and Mihalic also elaborated on the difficulty concerning a
conservative traction boundary condition when related to quantities defined in
current configuration. Applying this approach to a biventricular heart model,
it was shown that the inclusion of pre-stressing leads to significant strain
contributions in the myocardium Peirlinck et al. (2018). In the work of
Hirschvogel et al. (2017) the imaging-based patient-specific ventricular model
was pre-stressed to low end-diastolic pressure to account for the imaged,
stressed configuration, using a modified updated Lagrangian FE formulation
proposed by Gee et al. (2010) where on element level only the deformation
gradient and the isoparametric Jacobian matrix was updated.

The modelling of the deformation response of biological tissue to the applied
loading additionally requires the determination of the parameters of the
nonlinear elastic material relations. Initially, patient-specific mechanical
properties of blood arteries were estimated by analytical equations such as
the Young-Laplace Law. However, the underlying assumptions of material
homogeneity and isotropy contradicts the known complexity of myocardial and
vascular tissue Moulton et al. (1995). Klotz et al. (2006) uncovered from ex
vivo and in vivo measurements of end diastolic pressure volume relations
(EDPVR) and subsequent regression analysis that normalized EDPVR curves of the
left ventricle (LV) had a consistent profile across humans and a variety of
mammals, regardless of etiology. Using the least-squared error between the so-
called ”Klotz curve” and the modelled ventricular pressure-volume curve as an
optimization objective increases the number of comparison points and
facilitates the fitting of non-linear myocardial material models, e.g. Genet
et al. (2014); Palit et al. (2018); Sack et al. (2018). The use of the ”Klotz
curve” requires at least on measured pressure-volume pair which is usually the
end diastolic pressure and volume. However, in cardiac mechanics specifically,
the in vivo patient-specific LV pressure is often not known as well, because
the only accurate method for LV pressure measurement is invasive
catheterization which is therefore not readily available on an individual
basis. Without prior knowledge of the so-called end-diastolic pressure-volume
relation, it is effectively not possible to calibrate the anisotropic highly
non-linear elastic material relations for myocardial tissue.

To address the problem when direct pressure measurements are unavailable, this
work proposes to combine a direct inverse method to determine the unloaded
configurations of biventricular patient-specific heart models with the
previously mentioned ”Klotz curve” and the bounded Levenberg-Marquart
parameter optimisation method . This allows for the simultaneous computation
of the unknown end-diastolic ventricular pressures, pre-stressing and elastic
material constants. In contrast to the related inverse unloaded configuration
determination approaches by Govindjee and Mihalic (1996) and Gee et al.
(2010), the variational principle in the known deformed configuration is
derived from the conventional Lagrangian formulation.

The paper is organized as follows: In Sec. 2 the basics of a standard total
Lagrangian approach are revised. Based on this, a non-standard total
Lagrangian approach is introduced in Sec. 3 which allows for direct
computation of the unloaded configuration based on a know loaded
configuration. Finally, in Sec. 5 the proposed method is applied to three
numerical examples modelling non-linear hyperelastic material.

##  2 Standard total Lagrangian approach

This section provides only in brief the basic principles of a total Lagrangian
approach including necessary kinematical assumptions, variational principle
and its governing equations. For further details the reader is referred to
e.g. Gurtin2010.

###  2.1 Kinematics

Let $\mathbb{E}(3)$ be the Euclidian vector space and
$\mathcal{B}\subset\mathbb{E}(3)$, where $\mathcal{B}$ is a three-dimensional
manifold defining a material body. A motion of $\mathcal{B}$ is represented by
a one parameter non-linear deformation mapping
$\boldsymbol{\varphi}_{t}:\mathcal{B}\rightarrow\mathcal{B}_{t}$, where
$t\in\mathbb{R}$ is the time and $\mathcal{B}_{t}$ is the current
configuration at time $t$. Accordingly, each material point
$\mathbf{X}\in\mathcal{B}$ is related to its placement $\mathbf{x}$ in the
spatial configuration $\mathcal{B}_{t}$ by the mapping

$\displaystyle\mathbf{x}=\boldsymbol{\varphi}\left(\mathbf{X},t\right)$ (1)

In what follows and without loss of generality we identify the body
$\mathcal{B}$ with its undeformed reference configuration at a fixed time
$t_{0}$. The deformation map possesses an invertible linear tangent map
$\mathbf{F}=\mbox{Grad}\,\boldsymbol{\varphi}$ denoted by the deformation
gradient, where the Jacobian $J=\mbox{det}\,\mathbf{F}>0$. The operator Grad
represents the gradient with respect to the reference configuration

$\displaystyle\mbox{Grad}:=\frac{\partial}{\partial\mathbf{X}}\,.$ (2)

The body $\mathcal{B}$ is parameterized by the Cartesian coordinates
$X_{i},\,i=1,\,2,\,3$. Here, and in what follows, Latin indices take the
values $1,\,2\,\,\mbox{or}\,\,3$. As a deformation measure we make use of the
right Cauchy-Green deformation tensor $\bf C$ defined by

$\displaystyle\mathbf{C}=\mathbf{F}^{T}\mathbf{F}\,.$ (3)

In the following we want to confine ourselves to the quasi static case.

###  2.2 Variational principle

Let us consider a non-linear boundary value problem on domain $\mathcal{B}$
with boundary $\partial\mathcal{B}$. Dirichlet boundary conditions are
prescribed on $\partial\mathcal{B}_{D}\subset\partial\mathcal{B}$ and Neumann
boundary conditions are prescribed on
$\partial\mathcal{B}_{N}=\partial\mathcal{B}\setminus\partial\mathcal{B}_{D}$.
We define the external virtual work in the Lagrangian form as follows

$\displaystyle\mathcal{W}_{ext}=\int_{\mathcal{B}}\mathbf{b}\cdot\delta\mathbf{u}\,\,dV+\int_{\partial\mathcal{B}_{N}}\mathbf{t}^{(\mathbf{n})}\cdot\delta\mathbf{u}\,\,dA$
(4)

where $\mathbf{b}$ denotes the body force, $\mathbf{t}^{(\mathbf{n})}$ the
external traction and $\mathbf{n}$ the outward normal on
$\partial\mathcal{B}_{N}$.

Furthermore, we assume now that the body under consideration $\mathcal{B}$ is
hyperelastic and possesses an elastic potential $\Psi$ represented by the
stored strain energy per unit mass $\psi(\mathbf{C})$. The variation of the
internal potential with respect to $\mathbf{C}$ in the Lagrangian form reads
as follows

$\displaystyle\delta\Psi=\int_{\mathcal{B}}\,\rho_{0}\,\frac{\partial\psi}{\partial\mathbf{C}}:\delta\mathbf{C}\,\,dV\,.$
(5)

In the static case and by considering only mechanical processes the first law
of thermodynamics provides the following variational statement

$\displaystyle\delta\Psi-\mathcal{W}_{ext}=0\,.$ (6)

Upon substitution of Eqs. (4) and (5) this variational statement can be
written as

$\displaystyle\int_{\mathcal{B}}\frac{1}{2}\,\mathbf{S}:\delta\mathbf{C}\,\,dV-\int_{\mathcal{B}}\mathbf{b}\cdot\delta\mathbf{u}\,\,dV-\int_{\partial\mathcal{B}_{N}}\,\mathbf{t}\cdot\delta\mathbf{u}\,\,dA=0\,,$
(7)

with the second Piola-Kirchhoff stress tensor $\mathbf{S}$ given by

$\displaystyle\mathbf{S}=2\rho_{0}\frac{\partial\psi}{\partial\mathbf{C}}\,.$
(8)

The governing equations of the above variational formulation (Eq. 7) can be
derived by applying Gauss’s divergence theorem and, assuming the variation
$\delta\mathbf{u}=\mathbf{0}$ on $\partial\mathcal{B}_{D}$ but otherwise
arbitrary, we have

$\displaystyle\mbox{Div}\left(\mathbf{F}\mathbf{S}\right)+\mathbf{b}=\mathbf{0}\quad\mbox{on}\,\,\mathcal{B}\,,$
(9)

and the corresponding Neumann boundary condition is given by

$\displaystyle\mathbf{F}\mathbf{S}\mathbf{n}=\mathbf{t}^{\left(\mathbf{n}\right)}\quad\mbox{on}\,\,\partial\mathcal{B}_{N}\,.$
(10)

where Div denotes the divergence operator in the reference configuration.
These field equations are supplemented by essential boundary conditions, the
so-called Dirichlet boundary conditions

$\displaystyle\mathbf{u}=\mathbf{g}\quad\mbox{on}\,\,\partial\mathcal{B}_{D}\,.$
(11)

with $\mathbf{g}$ being the displacement field prescribed on
$\partial\mathcal{B}_{D}$.

##  3 Non-standard total Lagrangian approach

If the unloaded and stress-free configuration $\mathcal{B}$ is not known but
only a certain deformed configuration $\hat{\mathcal{B}}$, the standard total
Lagrangian approach as introduced in the previous section cannot be directly
used but needs to be modified. The reason for this is found in the fact that
in the total Lagrangian approach the displacement field and its gradient as
well as stress and strain quantities refer to the reference configuration
which is unknown.

Figure 1: Configuration spaces.

In order to proceed let us first introduce some deformation maps: As proposed
by Gee et al. (2010) the unknown reference configuration $\mathcal{B}$ is
linked with the know deformed configuration $\hat{\mathcal{B}}$ by
$\hat{\mathbf{F}}$ and any subsequent deformed configuration $\mathcal{B}_{t}$
by $\bar{\mathbf{F}}$. Thus,

$\displaystyle\mathbf{F}=\bar{\mathbf{F}}\,\hat{\mathbf{F}}=\frac{\partial\mathbf{x}}{\partial\hat{x}_{i}}\frac{\partial\hat{x}_{i}}{\partial\mathbf{X}}$
(12)

relates $\mathcal{B}$ with $\mathcal{B}_{t}$ as illustrated in Fig. 1.

###  3.1 Inverse problem

First, the objective is to find

$\displaystyle\hat{\mathbf{F}}=\frac{\partial\hat{\mathbf{x}}}{\partial\mathbf{X}}$
(13)

with $\hat{\mathbf{x}}\in\hat{\mathcal{B}}$ being known which we want to call
”inverse problem”. For this the standard variational principle needs to be
reformulated such as to determine $\mathbf{X}\in\mathcal{B}$ by employing the
inverse of the tangent map (Eq. 13)

$\displaystyle\hat{\mathbf{F}}^{-1}=\frac{\partial(\hat{\mathbf{x}}-\hat{\mathbf{u}})}{\partial\hat{\mathbf{x}}}\,,$
(14)

where the displacement field $\hat{\mathbf{u}}$ translates $\mathcal{B}$ to
the $\hat{\mathcal{B}}$.

Let us start with the standard variational formulation in its Lagrangian form
(Eq. 7) evaluated when $\mathcal{B}$ is deformed into $\hat{\mathcal{B}}$ such
that it holds $\mathbf{F}=\hat{\mathbf{F}}$:

$\displaystyle\int_{\mathcal{B}}\,\hat{\mathbf{F}}\mathbf{S}:\delta\hat{\mathbf{F}}\,dV-\mathcal{W}_{ext}=0$
(15)

where the second Piola-Kirchhoff stress tensor is now a function of
$\hat{\mathbf{C}}=\hat{\mathbf{F}}^{T}\hat{\mathbf{F}}$, i.e.

$\displaystyle\mathbf{S}=2\rho_{0}\frac{\partial\psi(\hat{\mathbf{C}})}{\partial\hat{\mathbf{C}}}\,.$
(16)

From

$\displaystyle\hat{\mathbf{F}}\hat{\mathbf{F}}^{-1}=\mathbf{1}$

which is varied with respect to $\hat{\mathbf{u}}$

$\displaystyle\delta\hat{\mathbf{F}}\hat{\mathbf{F}}^{-1}+\hat{\mathbf{F}}\delta\hat{\mathbf{F}}^{-1}=\mathbf{0}$

we obtain

$\displaystyle\delta\hat{\mathbf{F}}=-\hat{\mathbf{F}}\delta\hat{\mathbf{F}}^{-1}\hat{\mathbf{F}}\,.$
(17)

This is substituted together with $dV=\hat{J}^{-1}d\hat{v}$ into the above
standard variational formulation (Eq. 15) and we have

$\displaystyle\int_{\hat{\mathcal{B}}}\,-\hat{J}^{-1}\hat{\mathbf{C}}\hat{\mathbf{S}}\hat{\mathbf{F}}^{T}:\delta\hat{\mathbf{F}}^{-1}\,d\hat{v}-\mathcal{W}_{ext}=0\,,$
(18)

which is with Eq. (14) and the Cauchy stress tensor given by
$\hat{\boldsymbol{\sigma}}=-\hat{J}^{-1}\hat{\mathbf{C}}\hat{\mathbf{S}}\hat{\mathbf{F}}^{T}$
analog to Gee et al. (2010); Govindjee and Mihalic (1996); Peirlinck et al.
(2018). Here, it is important to note that, as the integration needs to be
performed over $\hat{\mathcal{B}}$, the approximation of geometry refers to
$\hat{\mathbf{x}}\in\hat{\mathcal{B}}$ and corresponding spatial derivatives
$\displaystyle\frac{\partial}{\partial\hat{\mathbf{x}}}\left(\bullet\right)$.

In terms of loading, we need discriminate between conservative and non-
conservative loading, i.e. whether it depends on deformation using the
standard total Lagriangian approach (Eq. 7). For non-conservative loading,
e.g. pressure loading $p$, using Nanson’s formula
$\hat{\mathbf{n}}\,d\hat{a}=\hat{J}\hat{\mathbf{F}}^{-T}\mathbf{n}\,dA$ which
establishes a relation between known deformed and corresponding unknown
undeformed surface elements $d\hat{\mathbf{a}}$ and $d\mathbf{A}$,
respectively, the external virtual work is expressed in the unknown undeformed
configuration as

$\displaystyle\mathcal{W}_{ext}=\int_{\partial\mathcal{B}_{N}}\,p\,\hat{J}\hat{\mathbf{F}}^{-T}\mathbf{n}\cdot\delta\hat{\mathbf{u}}\,dA$
(19)

whereas in the known deformed configuration as

$\displaystyle\mathcal{W}_{ext}=\int_{\partial\hat{\mathcal{B}}_{N}}\,p\,\hat{\mathbf{n}}\cdot\delta\hat{\mathbf{u}}\,d\hat{a}$
(20)

with $\hat{\mathbf{n}}$ denoting the normal vector on
$\partial\hat{\mathcal{B}}$ which does not depend on the deformation. The
latter is a crucial point to realize, as it is in contrast to standard total
Lagriangian approach. When dealing with the inverse problem
$\hat{\mathcal{B}}$ is assumed fixed whereas $\mathcal{B}_{0}$ needs to be
incrementally approached due to the nonlinear nature of the problem.
Consequently, when solving the the variational problem (Eq. 18) numerically,
e.g. employing the Newton-Raphson method, the pressure does not appear in the
variational statement’s linearized form. Conversely, when applying a
conservative loading, e.g. a traction force acting on the undeformed surface
$\partial\mathcal{B}_{N}$, we have

$\displaystyle\mathcal{W}_{ext}=\int_{\partial\mathcal{B}_{N}}\,\mathbf{t}\cdot\delta\hat{\mathbf{u}}\,dA$
(21)

or with respect to the known deformed surface $\partial\hat{\mathcal{B}}_{N}$

$\displaystyle\int_{\partial\hat{\mathcal{B}}_{N}}\,\mathbf{t}\cdot\delta\hat{\mathbf{u}}\,\hat{J}^{-1}\,|\hat{\mathbf{F}}^{T}\,\hat{\mathbf{n}}|\,d\hat{a}$
(22)

where the sizes of undeformed and deformed surface elements can be obtained
form Nanson’s formula
$|\mathbf{n}|\,dA=\hat{J}^{-1}\,|\hat{\mathbf{F}}^{T}\,\hat{\mathbf{n}}|\,d\hat{a}$.
Consequently, for inverse problem the linearization of the surface traction
contribution (Eq. 22) does not vanish.

###  3.2 Forward problem

Once the map $\hat{\mathbf{F}}$ linking unknown undeformed configuration
$\mathcal{B}_{0}$ and known deformed configuration $\hat{\mathcal{B}}$ has
been determined as described in Sec. 3.1, any subsequent loading application
which we want to term ”forward problem” can be dealt with computing the total
deformation gradient according to
$\mathbf{F}=\bar{\mathbf{F}}\hat{\mathbf{F}}$ where $\bar{\mathbf{F}}$ is a
function of the displacement increment $\bar{\mathbf{u}}$ due to some
additional loading and $\hat{\mathbf{F}}$ remains constant. In practise, when
employed to a numerical method, $\hat{\mathbf{F}}$ needs to be stored at all
points of interest, e.g. numerical integration points.

Now, for forward problem the original variational formulation (Eq. 7) is
reformulated with $dV=\hat{J}^{-1}d\hat{v}$ as follows:

$\displaystyle\int_{\hat{\mathcal{B}}}\,\hat{J}^{-1}\mathbf{F}\mathbf{S}\hat{\mathbf{F}}^{T}:\delta\bar{\mathbf{F}}\,d\hat{v}-\mathcal{W}_{ext}=0$
(23)

where the variation of $\mathbf{F}$ expressed as

$\displaystyle\delta\mathbf{F}=\delta\bar{\mathbf{F}}\hat{\mathbf{F}}\,.$ (24)

For the latter it is important to realize that, as $\hat{\mathbf{F}}$ is
constant during the subsequent ”forward” computation, the variation is only
performed with respect $\bar{\mathbf{u}}=\mathbf{x}-\hat{\mathbf{x}}$, as
$\bar{\mathbf{F}}=\bar{\mathbf{F}}(\bar{\mathbf{u}})$ is linking
$\hat{\mathcal{B}}$ and $\mathcal{B}_{t}$.

With regards to deformation dependent loading we consider

$\displaystyle\mathcal{W}_{ext}=\int_{\partial\hat{\mathcal{B}}_{N}}\,p\,\bar{J}\bar{\mathbf{F}}^{-T}\hat{\mathbf{n}}\cdot\delta\bar{\mathbf{u}}\,d\hat{a}$
(25)

using Nanson’s formula
$\bar{\mathbf{n}}\,d\bar{a}=\bar{J}\bar{\mathbf{F}}^{-T}\hat{\mathbf{n}}\,d\hat{a}$
which establishes the relation between unknown deformed and corresponding
known deformed surface elements $d\bar{\mathbf{a}}$ and $d\hat{\mathbf{a}}$,
respectively. For deformation independent traction loading (Eq. (22)) still
applies. Consequently, for the forward problem, as for the standard total
Lagrangian approach, the pressure contribution (Eq. (25)) is to be linearized
with respect to $\bar{\mathbf{u}}$ whereas the linearization of the surface
traction contribution does vanish, as $\hat{\mathbf{F}}$ is now constant.

##  4 Inverse parameter optimisation method

As previously mentioned, the determination of elastic myocardial material
parameters requires knowledge of ventricular filling pressure and
corresponding deformation. The latter can include deformation metrics such as
change of cavity volume, short-axis and long-axis cavity expansion, twisting
about the long-axis or globally and regionally averaged myocardial strains,
e.g. Genet et al. (2014); Palit et al. (2018); Sack et al. (2018).

Here, we aim at determining the elastic myocardial material parameters of
biventricular (BV) models without prior knowledge of ventricular filling
pressure data by exploiting the characteristic exponential nature of the
”Klotz curve” fitted to three patient-specific left ventricular volume data,
namely ESV and EDV from Cardiovascular Magnetic Resonance (CMR) imaging as
well as the undeformed cavity volume, $V_{0}$, as obtained by inverse
computation using the non-standard total Lagrangian approach presented in Sec.
3 implemented in a finite element method (FEM)-based code.

This is achieved by iterative parameter optimisation employing the Levenberg-
Marquardt algorithm (LVM) Guyon and Le Riche (2000). Further details regarding
the specific LVM implementation can be found in Rama and Skatulla (2019). This
iterative scheme is initialized by assuming a physiological EDP and elastic
myocardial material parameters for the LV which serves to construct an initial
volume-normalized ”Klotz curve” Klotz et al. (2006). Before elaborating on the
details of the proposed iterative scheme, the important relations needed to
set-up the ”Klotz curve” are briefly summarized. Given a measured pair of LV
cavity volume, $V_{m}$, and corresponding cavity pressure, $P_{m}$, the
unloaded LV cavity volume can be estimated as

$V_{0}=V_{m}\,\left(0.6-0.006\,P_{m}\right)$ (26)

and at 30 mmHg as

$V_{30}=V_{0}+\frac{V_{m,n}-V_{0}}{(P_{m}/A_{n})^{(1/B_{n})}}$ (27)

where the normalized measured volume is given by

$V_{m,n}=\frac{V_{m}-V_{0}}{V_{30}-V_{0}}$ (28)

and for humans the coefficients $A_{n}=27.8$ and $B_{n}=2.76$. Then the EDPVR
is formulated as

$P_{LV}=\alpha\,V_{LV}^{\beta}$ (29)

with

$\alpha=\displaystyle\frac{30}{\left(V_{30}\right)^{\text{Log}(P_{m}/30)/\text{Log}(V_{m}/V_{30})}}$
(30)

and

$\beta=\displaystyle\frac{\text{Log}(P_{m}/30)}{\text{Log}(V_{m}/V_{30})}\,.$
(31)

Now, in the absence of LV pressure measurements, EDP is to be iteratively
estimated taking $P_{m}=\text{EDP}$ and $V_{m}=\text{EDV}$. For the initial
setup of the ”Klotz curve”, a physiologically reasonable EDP value needs to be
chosen. ESP and $V_{0}$ can then be readily obtained via Eq. (29) knowing ESV
from CMR and $P_{0}=0$. The LVM cost function considers the fitting error of
$V_{0}$ and EDV between ”Klotz curve” and the FEM solution. As cavity pressure
data are generally not available for the RV, it is therefore assumed that the
”Klotz curve” (Eq. (29)) can provide an indication of the EDPVR of the RV as
well. Due to this limitation, the LVM cost function only includes the error
regarding EDV but not $V_{0}$ for the RV. For each subsequent LVM iteration,
updated elastic parameter and EDP values are obtained for LV and RV from the
optimisation algorithm allowing for the update of the ”Klotz curve” and
subsequently ESP and $V_{0}$ values as well. A noteworthy aspect of this
approach is that the fitting target, which is the ”Klotz curve”, is
continuously changing which is not the case in the standard application of
LVM. The iterative parameter optimisation is terminated once the cost function
error falls below a given threshold and the root-square error of the entire
EDPVR of ”Klotz curve” and FE simulation is smaller than 2.00 mmHg.

Figure 2: Inverse procedure combining ”Klotz Curve”, FE simulation and
Levenberg-Marquardt (LVM) optimization schemes.

##  5 Numerical examples

The following examples presented in this section aim to demonstrate the
potential of the proposed inverse parameter optimisation method outlined in
Sec. 4. In this study, only the diastolic filling of the biventricular (BV)
cavities will be examined.

Table 1: Ventricular cavity volumes measured by CMR. | Left ventricle | Right ventricle
---|---|---
| ESV | EDV | SV | ESV | EDV | SV
Case 1 | 90.00 | 168.00 | 78.00 | 61.00 | 96.00 | 78.00
Case 2 | 86.00 | 170.00 | 84.00 | 67.00 | 118.00 | 51.00
Case 3 | 82.00 | 157.00 | 75.00 | 78.00 | 160.00 | 82.00

We constructed patient-specific heart models of three healthy patients that
were scanned using Cardiovascular Magnetic Resonance (CMR) imaging at the
Groote Schuur Hospital, University of Cape Town, South Africa. University of
Cape Town, Faculty of Health Sciences Human Research Ethics Committee approval
(REF: 686/2018) and patients’ consent were obtained to conduct research on
unidentified human data. The left and right ventricular cavity volumes at end-
systole and end-diastole, respectively, are presented in Tab. 1. The
anatomical models were meshed using 3326, 2765 and 3246 linear tetrahedral
elements, for the cases 1, 2, and 3, respectively. ESP cavity filling pressure
boundary conditions were applied to the endocardial surfaces of LV and RV
according to Eq. (20) computing first the unloaded configuration as outlined
in Sec. 3.1 and EDP cavity filling pressure boundary conditions according to
Eq. (25) computing subsequently the end diastolic configuration as outlined in
Sec. 3.2. In terms of the Dirichlet boundary conditions, the surface of the
heart’s base is fixed to restrain vertical direction movement. To allow for
torsional behaviour and wall thickening and thinning, a Dirichlet boundary
condition is weakly imposed through application of an elastic line force of
$0.1kN$ acting in tangential direction around the epicardial base. The elastic
line forces effectively prevent rigid body motion in short-axis directions as
well.

The material behaviour of the BV heart muscle tissue is assumed to be
nonlinear, anisotropic and nearly incompressible described by the following
orthotropic strain energy function Usyk et al. (2000)

$\displaystyle\psi$ $\displaystyle=$
$\displaystyle\frac{A}{2}\left(\exp^{BQ}-1\right)+A_{\text{comp}}\Big{(}\text{det}J\;\text{ln}\left(\text{det}J\right)-\text{det}J+1\Big{)},$
(32)

where parameters $A$ and $B$ are stiffness factors and
$J=\text{det}\mathbf{F}$ is the Jacobian quantifying volumetric deformation of
cardiac tissue is linked to the constant $A_{comp}$ controlling the
compressibility of the myocardium. Additionally, the exponent $Q$ is defined
in terms of the Green strain tensor and the material directions defining
structural tensors $\mathbf{M}_{f}$, $\mathbf{M}_{s}$ and $\mathbf{M}_{n}$ as
follows Rama et al. (2016):

$\displaystyle Q$ $\displaystyle:=$ $\displaystyle
a_{1}(\text{tr}(\mathbf{M}_{f}\mathbf{E}))^{2}+a_{2}(\text{tr}(\mathbf{M}_{s}\mathbf{E}))^{2}+a_{3}(\text{tr}(\mathbf{M}_{n}\mathbf{E}))^{2}$
(33)
$\displaystyle+a_{4}\,\text{tr}(\mathbf{M}_{f}\mathbf{E}^{2})+a_{5}\,\text{tr}(\mathbf{M}_{s}\mathbf{E}^{2})+a_{6}\,\text{tr}(\mathbf{M}_{n}\mathbf{E}^{2}),$

where $a_{i},\,i=1,...,6$ are the anisotropy coefficients associated with the
three preferred material directions, namely fiber axis, $\mathbf{V}_{f}$,
sheet axis, $\mathbf{V}_{s}$, and sheet normal axis, $\mathbf{V}_{n}$. These
vectors construct an orthonormal basis and allow for the formulation of the
so-called structural tensors

$\displaystyle\mathbf{M}_{f}=\mathbf{V}_{f}\otimes\mathbf{V}_{f},\mathbf{\ }\
M_{s}=\mathbf{V}_{s}\otimes\mathbf{V}_{s},\ \
\mathbf{M}_{n}=\mathbf{V}_{n}\otimes\mathbf{V}_{n}.$ (34)

The values of the anisotropy material constants used are listed in Table 2.

Table 2: Anisotropy coefficients used. $A_{comp}$ (kPa) | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | $a_{6}$
---|---|---|---|---|---|---
100 | 6 | 7 | 3 | 12 | 3 | 3

The varying fibre directions throughout the LV and RV walls are computed
adopting an algorithm developed by Wong and Kuhl Wenk et al. (2011) but using
moving least square (MLS) based approximations as proposed by Skatulla and
Sansour (2016) instead of the solving a Poisson problem. The resulting fibre
direction distribution is shown for Case 2 in Fig. 3. The prescribed fibre
direction angles at the endocardium and epicardium of both LV and RV are
$+72^{o}$ and $-57^{o}$, respectively Rama and Skatulla (2019).

Figure 3: Myocardial fibre orientation of a bi-ventricular model.

###  5.1 Iterative parameter estimation

The iterative parameter estimation approach is used to obtain patient-specific
values for the characteristic cavity pressures $EDP$ and $ESP$, respectively,
as well as the stiffness parameters $A$ and $B$ for LV and RV. The
corresponding physiologically motivated initial values are listed in Tab. 3.
The LVM scaling parameters are set as $\lambda^{\text{LVM}}=1$ and
$\nu^{\text{LVM}}=2$, respectively and the stopping criterion of the cost
function residual is chosen as $\gamma^{\text{LVM}}_{\epsilon}=1X10^{-2}$.

Table 3: Initial material parameter and cavities pressure values in kPa Zhou
et al. (2020); Mielniczuk et al. (2007); Wang et al. (2013)

.  Left ventricle Right ventricle A (kPa) B $EDP$ A (kPa) B $EDP$ Case 1 0.10
1.00 2.00 0.10 1.00 0.50 Case 2 0.10 1.00 3.00 0.10 1.00 1.50 Case 3 0.10 1.00
2.00 0.10 1.00 2.50

The final obtained stiffness parameter and characteristic pressure values are
listed in Tab. 4.

Table 4: Final estimated material parameters and cavities pressure values in kPa. | Left ventricle | Right ventricle
---|---|---
| A (kPa) | B | $ESP$ | $EDP$ | A (kPa) | B | $ESP$ | $EDP$
Case 1 | 0.10 | 1.25 | 0.04 | 1.78 | 0.11 | 1.05 | 0.08 | 0.85
Case 2 | 0.09 | 1.20 | 0.04 | 2.55 | 0.12 | 1.21 | 0.04 | 1.06
Case 3 | 0.09 | 1.07 | 0.04 | 1.72 | 0.10 | 0.88 | 0.02 | 1.25

The converged iterative parameter optimisation scheme gave for the unloaded LV
cavity volumes $V_{0}=86.66,\,83.32\,\text{and}\,81.94\,\text{mL}$ and for the
unloaded RV cavity volumes $V_{0}=53.92,\,65.18\,\text{and}\,86.99\,\text{mL}$
for cases 1, 2 and 3, respectively. Figs. 4 \- 6 show the fitted EDPVR curves
illustrating the FEM solutions by red lines and the targeted ”Klotz curve” by
blue circles indicating good agreement. The root square mean error (RMSE) are
0.225, 0.225 and 0.375 mmHg for cases 1, 2 and 3 respectively at the LV and
0.150, 0.150 and 0.075 mmHg also for cases 1, 2 and 3 respectively.

Figure 4: EDPVR curve predicted by Usyk model (red) and method of Klotz et al.
(blue) as a physiological benchmark for case 1. Figure 5: EDPVR curve
predicted by Usyk model (red) and method of Klotz et al. (blue) as a
physiological benchmark for case 2. Figure 6: EDPVR curve predicted by Usyk
model (red) and method of Klotz et al. (blue) as a physiological benchmark for
case 3.

The average percentage difference between the simulated and targeted CMR
LVEDVs was approximately 0.0007$\%$, and 0.0020$\%$ for the RVEDVs. Similarly,
the average percentage difference between the simulated and targeted Klotz
$V_{0}$ for the LV was approximately 0.00005$\%$ and 5.3314$\%$ for the RV.
The observed higher percentage error in the RV initial volume confirm that we
did not target RV initial volume in our calibration. For the segmented and CMR
ESV, there is 0.92$\%$ average error difference at the LV and 3.00$\%$ at the
RV, detailed result is on Table. 5.

Table 5: Percentage error of LV and RV volumetric targets and results from the calibration. | Left ventricle | Right ventricle
---|---|---
| $V_{0}$ | LVEDV | LVESV | $V_{0}$ | RVEDV | RVESV
Case 1 | 0.00005 | 0.00172 | 1.42 | 1.73 | 0.0003 | 3.92
Case 2 | 0.00006 | 0.00029 | 0.01 | 0.01 | 0.0020 | 1.86
Case 3 | 0.00002 | 0.00005 | 1.33 | 14.25 | 0.0036 | 3.22

###  5.2 Stress and strain distributions

The contour plots in Figs. 7 and 8 reveal high myofibre stress and strain
concentrations at the base which is to be expect due to the geometrical
support. The global residual fibre stress, separately averaged over LV and RV
myocardium, at end systole is found as 0.013, 0.018 and 0.006 kPa for LV at
the longitudinal, circumferential and radial respectively and 0.031, 0.038,
and 0.012 kPa for RV longitudinal, circumferential and radial respectively. A
detailed analysis of regional strain in the LV with respect to the global
longitudinal, circumferential and radial directions is presented. Endocardial
strain comparison demonstrated moderate agreement of global strains,
separately averaged over the LV and RV myocardium, respectively. Global
longitudinal strain values (GLS) for cases 1 to 3 are 5.75, 6.90 and 7.14$\%$
for the FE model simulation, while the CMR in vivo recorded values for cases 1
to 3 are given as 12.30, 17.80 and 21.60$\%$ respectively. The global
circumferential strain values (GCS) for the FE simulations are 22.2, 25.82 and
25.82$\%$, and 20.20, 20.40 and 25.70$\%$ for the in vivo. Similarly, the
global radial strain values (GRS) for the FE simulations are -15.84, -17.37,
and -16.39$\%$, and -13.90, -14.00, and -16.20$\%$ in vivo for the cases 1, 2
and 3 respectively. The GCS and GRS are well comparable in both FE model
simulation and the CMR in vivo, but this was not the case with longitudinal
strains, whereby a regional analysis of the longitudinal strains revealed poor
agreement between simulated and measured in vivo strains. Although
discrepancies have been noted between the FEM and in vivo longitudinal strain
Sack et al. (2018) due to the absence of data on the true distribution of the
patient-specific myofiber orientation angles, which affects the
circumferential-longitudinal compliance ratio. Also the CMR long-axis scan
providing longitudinal strain and the short-axis scan providing
circumferential and radial strains are obtained from different cardiac cycles
whereas the FEM calculate these strains using the same cardiac cycle.

Figure 7: Fibre stress at the end-diastole. Figure 8: Fibre strain at the end-
diastole.

##  6 Conclusion

In this paper, a direct inverse method to determine the unloaded
configurations of biventricular patient-specific heart models is combined with
the bounded Levenberg-Marquart parameter optimisation method to estimated
besides elastic material parameters also the unknown end-systolic and end-
diastolic ventricular pressures in silico. The unloaded configuration is
computed by a modified total Lagrangian approach as implemented in standard
FEM which implicitly provides for the incorporation of pres-stressing and pre-
straining of the myocardium as linked to the end-systolic pressure. The
proposed variational principle is derived from the standard Lagrangian
formulation, in contrast to the similar inverse unloaded configuration
determination methods by Govindjee and Mihalic (1996) and Gee et al. (2010).
The unknown characteristic ventricular pressures and elastic material
parameters can then be iteratively computed by taking advantage of the highly
nonlinear exponential nature of the EDPVR and fitting the patient-specific LV
EDPVR to the analytical volume-normalized ”Klotz EDPVR curve”. Making use of
CMR-derived anatomical heart models, the iterative framework therefore
provides the means to non-invasively estimate LV pressure, pre-strain and
elastic material parameters exploiting their biomechanical interdependence.

Our results indicate good agreement between the EDPVR curves obtained from FE
simulation and given by the ”Klotz Curve”. Additionally, there is a close
agreement in the FE global circumferential and radial strains when compared
with the CMR global strain data. The global longitudinal strain shows poor
agreement between FE and CMR measured strains. This discrepancy has been
attributed to the absence of data on the true distribution of the patient-
specific myofiber orientation distribution, which affects the circumferential-
longitudinal compliance ratio Sack et al. (2018). Also the CMR long-axis scan
providing longitudinal strain and short-axis scan providing circumferential
and radial strains are obtained from different cardiac cycles whereas FEM
calculates these strains within the same cardiac cycle. Another limitation of
this study is the absence of in vivo pressure measurements to validate our
results. In future work, we shall obtain pressure data in vivo to perform
validation of our proposed method, also patient-specific fiber orientation
distribution data from DTMRI will help to mitigate the error in FE and CMR
global longitudinal strains.

## References

## References

  * Baillargeon et al. [2014] Brian Baillargeon, Nuno Rebelo, David D Fox, Robert L Taylor, and Ellen Kuhl.  The living heart project: a robust and integrative simulator for human heart function.  _European Journal of Mechanics-A/Solids_ , 48:38–47, 2014\.
  * Bols et al. [2013] Joris Bols, Joris Degroote, Bram Trachet, Benedict Verhegghe, Patrick Segers, and Jan Vierendeels.  A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels.  _Journal of computational and Applied mathematics_ , 246:10–17, 2013.
  * Gee et al. [2010] Michael W Gee, Ch Förster, and WA Wall.  A computational strategy for prestressing patient-specific biomechanical problems under finite deformation.  _International Journal for Numerical Methods in Biomedical Engineering_ , 26(1):52–72, 2010.
  * Genet et al. [2014] Martin Genet, Lik Chuan Lee, Rebecca Nguyen, Henrik Haraldsson, Gabriel Acevedo-Bolton, Zhihong Zhang, Liang Ge, Karen Ordovas, Sebastian Kozerke, and Julius 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(2):142–152, 2014.
  * Govindjee and Mihalic [1996] Sanjay Govindjee and Paul A Mihalic.  Computational methods for inverse finite elastostatics.  _Computer Methods in Applied Mechanics and Engineering_ , 136(1-2):47–57, 1996.
  * Gray and Pathmanathan [2018] Richard A Gray and Pras Pathmanathan.  Patient-specific cardiovascular computational modeling: diversity of personalization and challenges.  _Journal of cardiovascular translational research_ , 11:80–88, 2018.
  * Guyon and Le Riche [2000] Frédéric Guyon and Rodolphe Le Riche.  Least squares parameter estimation and the levenberg-marquardt algorithm: Deterministic analysis, sensitivities and numerical experiments.  Technical Report 041/99, Institut National des Sciences Appliquées, 2000.
  * Hirschvogel et al. [2017] Marc Hirschvogel, Marina Bassilious, Lasse Jagschies, Stephen M Wildhirt, and Michael 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(8):e2842, 2017.
  * Klotz et al. [2006] Stefan Klotz, Ilan Hay, Marc L Dickstein, Geng-Hua Yi, Jie Wang, Mathew S Maurer, David A Kass, and Daniel Burkhoff.  Single-beat estimation of end-diastolic pressure-volume relationship: a novel method with potential for noninvasive application.  _American Journal of Physiology-Heart and Circulatory Physiology_ , 291(1):H403–H412, 2006.
  * Mielniczuk et al. [2007] Lisa M Mielniczuk, Gervasio A Lamas, Greg C Flaker, Gary Mitchell, Sidney C Smith, Bernard J Gersh, Scott D Solomon, Lemuel A Moyé, Jean L Rouleau, John D Rutherford, et al.  Left ventricular end-diastolic pressure and risk of subsequent heart failure in patients following an acute myocardial infarction.  _Congestive Heart Failure_ , 13(4):209–214, 2007\.
  * Moulton et al. [1995] Michael J Moulton, Lawrence L Creswell, Ricardo L Actis, Kent W Myers, Michael W Vannier, Barna A Szabo, and Michael K Pasque.  An inverse approach to determining myocardial material properties.  _Journal of biomechanics_ , 28(8):935–948, 1995\.
  * Nikou et al. [2016] Amir Nikou, Shauna M Dorsey, Jeremy R McGarvey, Joseph H Gorman, Jason A Burdick, James J Pilla, Robert C Gorman, and Jonathan F Wenk.  Computational modeling of healthy myocardium in diastole.  _Annals of biomedical engineering_ , 44:980–992, 2016.
  * Omens et al. [2003] JH Omens, AD McCulloch, and JC Criscione.  Complex distributions of residual stress and strain in the mouse left ventricle: experimental and theoretical models.  _Biomechanics and modeling in mechanobiology_ , 1(4):267–277, 2003.
  * Palit et al. [2018] Arnab Palit, Sunil K Bhudia, Theodoros N Arvanitis, Glen A Turley, and Mark A Williams.  In vivo estimation of passive biomechanical properties of human myocardium.  _Medical & biological engineering & computing_, 56:1615–1631, 2018.
  * Peirlinck et al. [2018] Mathias Peirlinck, Matthieu De Beule, Patrick Segers, and Nuno Rebelo.  A modular inverse elastostatics approach to resolve the pressure-induced stress state for in vivo imaging based cardiovascular modeling.  _Journal of the mechanical behavior of biomedical materials_ , 85:124–133, 2018.
  * Pourmodheji et al. [2021] Reza Pourmodheji, Zhenxiang Jiang, Christopher Tossas-Betancourt, C Alberto Figueroa, Seungik Baek, and Lik-Chuan Lee.  Inverse modeling framework for characterizing patient-specific microstructural changes in the pulmonary arteries.  _Journal of the mechanical behavior of biomedical materials_ , 119:104448, 2021.
  * Rajagopal et al. [2007] Vijayaraghavan Rajagopal, Jae-Hoon Chung, David Bullivant, Poul MF Nielsen, and Martyn P Nash.  Determining the finite elasticity reference state from a loaded configuration.  _International Journal for Numerical Methods in Engineering_ , 72(12):1434–1451, 2007.
  * Rama et al. [2016] Ritesh R Rama, Sebastian Skatulla, and Carlo Sansour.  Real-time modelling of diastolic filling of the heart using the proper orthogonal decomposition with interpolation.  _International Journal of Solids and Structures_ , 96:409–422, 2016.
  * Rama and Skatulla [2019] Ritesh Rao Rama and Sebastian Skatulla.  Real-time nonlinear solid mechanics computations for fast inverse material parameter optimization in cardiac mechanics.  _Journal of Engineering Mechanics_ , 145(4):04019020, 2019.
  * Sack et al. [2018] Kevin L Sack, Eric Aliotta, Daniel B Ennis, Jenny S Choy, Ghassan S Kassab, Julius M Guccione, and Thomas Franz.  Construction and validation of subject-specific biventricular finite-element models of healthy and failing swine hearts from high-resolution dt-mri.  _Frontiers in Physiology_ , 9:539, 2018.
  * Shi et al. [2019] Xiaodan Shi, Yue Liu, Katherine M Copeland, Sara R McMahan, Song Zhang, J Ryan Butler, Yi Hong, Michael Cho, Pietro Bajona, Huajian Gao, et al.  Epicardial prestrained confinement and residual stresses: a newly observed heart ventricle confinement interface.  _Journal of the Royal Society Interface_ , 16(152):20190028, 2019.
  * Skatulla and Sansour [2016] S Skatulla and C Sansour.  On a path-following method for non-linear solid mechanics with applications to structural and cardiac mechanics subject to arbitrary loading scenarios.  _International Journal of Solids and Structures_ , 96:181–191, 2016.
  * Usyk et al. [2000] TP Usyk, R Mazhari, and AD McCulloch.  Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle.  _Journal of elasticity and the physical science of solids_ , 61(1):143–164, 2000.
  * Wang et al. [2013] HM Wang, H Gao, XY Luo, C Berry, BE Griffith, RW Ogden, and TJ Wang.  Structure-based finite strain modelling of the human left ventricle in diastole.  _International journal for numerical methods in biomedical engineering_ , 29(1):83–103, 2013.
  * Wenk et al. [2011] Jonathan F Wenk, Parastou Eslami, Zhihong Zhang, Chun Xu, Ellen Kuhl, Joseph H Gorman III, J Daniel Robb, Mark B Ratcliffe, Robert C Gorman, and Julius M Guccione.  A novel method for quantifying the in-vivo mechanical effect of material injected into a myocardial infarction.  _The Annals of thoracic surgery_ , 92(3):935–941, 2011.
  * Zhou et al. [2020] Xia Zhou, Mingming Lei, Donghui Zhou, Guoqing Li, Zhiying Duan, Shaobo Zhou, and Yuanzhe Jin.  Clinical factors affecting left ventricular end-diastolic pressure in patients with acute st-segment elevation myocardial infarction.  _Ann Palliat Med_ , 9(4):1834–1840, 2020.
  * Zhuan and Luo [2022] X Zhuan and XY Luo.  Volumetric growth of soft tissues evaluated in the current configuration.  _Biomechanics and Modeling in Mechanobiology_ , 21(2):569–588, 2022.