A Self-Taught Artificial Agent for Multi-Physics
Computational Model Personalization
Dominik Neumanna,c,∗, Tommaso Mansib , Lucian Itud,e ,
Bogdan Georgescub , Elham Kayvanpourf , Farbod Sedaghat-Hamedanif ,
Ali Amrf , Jan Haasf , Hugo Katusf , Benjamin Mederf , Stefan Steidlc ,
Joachim Horneggerc , Dorin Comaniciub

arXiv:1605.00303v1 [] 1 May 2016

a

Medical Imaging Technologies, Siemens Healthcare GmbH, Erlangen, Germany
b
Medical Imaging Technologies, Siemens Healthcare, Princeton, USA
c
Pattern Recognition Lab, FAU Erlangen-Nürnberg, Erlangen, Germany
d
Siemens Corporate Technology, Siemens SRL, Brasov, Romania
e
Transilvania University of Brasov, Brasov, Romania
f
Department of Internal Medicine III, University Hospital Heidelberg, Germany

Abstract
Personalization is the process of fitting a model to patient data, a critical
step towards application of multi-physics computational models in clinical
practice. Designing robust personalization algorithms is often a tedious,
time-consuming, model- and data-specific process. We propose to use artificial intelligence concepts to learn this task, inspired by how human experts
manually perform it. The problem is reformulated in terms of reinforcement learning. In an off-line phase, Vito, our self-taught artificial agent,
learns a representative decision process model through exploration of the
computational model: it learns how the model behaves under change of
parameters. The agent then automatically learns an optimal strategy for
on-line personalization. The algorithm is model-independent; applying it
to a new model requires only adjusting few hyper-parameters of the agent
and defining the observations to match. The full knowledge of the model
itself is not required. Vito was tested in a synthetic scenario, showing that
it could learn how to optimize cost functions generically. Then Vito was
applied to the inverse problem of cardiac electrophysiology and the personalization of a whole-body circulation model. The obtained results suggested
that Vito could achieve equivalent, if not better goodness of fit than standard methods, while being more robust (up to 11% higher success rates)
and with faster (up to seven times) convergence rate. Our artificial intelligence approach could thus make personalization algorithms generalizable
and self-adaptable to any patient and any model.
∗

Corresponding author
Email address: dominik.neumann@siemens.com (Dominik Neumann)

Preprint submitted to Medical Image Analysis

May 3, 2016

Keywords: Computational Modeling, Model Personalization,
Reinforcement Learning, Artificial Intelligence.

1. Introduction
Computational modeling attracted significant attention in cardiac research over the last decades (Clayton et al., 2011; Frangi et al., 2001; Hunter
and Borg, 2003; Kerckhoffs et al., 2008; Krishnamurthy et al., 2013; Kuijpers et al., 2012; Noble, 2002). It is believed that computational models
can improve patient stratification and therapy planning. They could become the enabling tool for predicting disease course and therapy outcome,
ultimately leading to improved clinical management of patients suffering
from cardiomyopathies (Kayvanpour et al., 2015). A crucial prerequisite for
achieving these goals is precise model personalization: the computational
model under consideration needs to be fitted to each patient. However, the
high complexity of cardiac models and the often noisy and sparse clinical
data still hinder this task.
A wide variety of manual and (semi-)automatic model parameter estimation approaches have been explored, including Aguado-Sierra et al. (2010,
2011); Augenstein et al. (2005); Chabiniok et al. (2012); Delingette et al.
(2012); Itu et al. (2014); Konukoglu et al. (2011); Le Folgoc et al. (2013);
Marchesseau et al. (2013); Neumann et al. (2014a,b); Prakosa et al. (2013);
Schmid et al. (2006); Seegerer et al. (2015); Sermesant et al. (2009); Wallman et al. (2014); Wang et al. (2009); Wong et al. (2015); Xi et al. (2013);
Zettinig et al. (2014). Most methods aim to iteratively reduce the misfit
between model output and measurements using optimization algorithms,
for instance variational (Delingette et al., 2012) or filtering (Marchesseau
et al., 2013) approaches. Applied blindly, those techniques could easily fail
on unseen data, if not supervised, due to parameter ambiguity, data noise
and local minima (Konukoglu et al., 2011; Neumann et al., 2014a; Wallman
et al., 2014). Therefore, complex algorithms have been designed combining
cascades of optimizers in a very specific way to achieve high levels of robustness, even on larger populations, i.e. 10 or more patients (Kayvanpour et al.,
2015; Neumann et al., 2014b; Seegerer et al., 2015). However, those methods
are often designed from tedious, trial-and-error-driven manual tuning, they
are model-specific rather than generic, and their generalization to varying
data quality cannot be guaranteed. On the contrary, if the personalization
task is assigned to an experienced human, given enough time, he almost always succeeds in manually personalizing a model for any subject (although
solution uniqueness is not guaranteed, but this is inherent to the problem).
There are several potential reasons why a human expert is often superior to standard automatic methods in terms of personalization accuracy
and success rates. First, an expert is likely to have an intuition of the

2

model’s behavior from his prior knowledge of the physiology of the modeled
organ. Second, knowledge about model design and assumptions, and model
limitations and implementation details certainly provide useful hints on the
“mechanics” of the model. Third, past personalization of other datasets
allows the expert to build up experience. The combination of prior knowledge, intuition and experience enables to solve the personalization task more
effectively, even on unseen data.
Inspired by humans and contrary to previous works, we propose to address the personalization problem from an artificial intelligence (AI) perspective. In particular, we apply reinforcement learning (RL) methods (Sutton
and Barto, 1998) developed in the AI community to solve the parameter
estimation task for computational physiological models. With its roots in
control theory on the one hand, and neuroscience theories of learning on the
other hand, RL encompasses a set of approaches to make an artificial agent
learn from experience generated by interacting with its environment. Contrary to standard (supervised) machine learning (Bishop, 2006), where the
objective is to compute a direct mapping from input features to a classification label or regression output, RL aims to learn how to perform tasks. The
goal of RL is to compute an optimal problem-solving strategy (agent behavior), e.g. a strategy to play the game “tic-tac-toe” successfully. In the AI
field, such a behavior is often represented as a policy, a mapping from states,
describing the current “situation” the agent finds itself in (e.g. the current
locations of all “X” and “O” on the tic-tac-toe grid), to actions, which allow
the agent to interact with the environment (e.g. place “X” on an empty cell)
and thus influence that situation. The key underlying principle of RL is that
of reward (Kaelbling et al., 1996), which provides an objective means for the
agent to judge the outcome of its actions. In tic-tac-toe, the agent receives
a high, positive reward if the latest action led to a horizontal, vertical or
diagonal row full of “X” marks (winning), and a negative reward (punishment) if the latest action would allow the opponent to win in his next move.
Based on such rewards, the artificial agent learns an optimal winning policy
through trial-and-error interactions with the environment.
RL was first applied to game (e.g. Tesauro, 1994) or simple control
tasks. However, the past few years saw tremendous breakthroughs in RL for
more complex, real-world problems (e.g. Barreto et al., 2014; Kveton and
Theocharous, 2012; Nguyen-Tuong and Peters, 2011). Some noteworthy examples include Mülling et al. (2013), where the control entity of a robot
arm learned to select appropriate motor primitives to play table tennis, and
Mnih et al. (2015), where the authors combine RL with deep learning to
train an agent to play 49 Atari games, yielding better performance than an
expert in the majority of them.
Motivated by these recent successes and building on our previous work
(Neumann et al., 2015), we propose an RL-based personalization approach,
henceforth called Vito, with the goal of designing a framework that can, for
3

Off-line phase

Assimilate model
behavior

On-line phase

Learn optimal
personalization
strategy

Personalize model
for new patients

Figure 1: Overview of Vito: a self-taught artificial agent for computational model personalization, inspired by how human operators approach the personalization problem.

the first time to our knowledge, learn by itself how to estimate model parameters from clinical data while being model-independent. As illustrated in
Fig. 1, first, like a human expert, Vito assimilates the behavior of the physiological model under consideration in an off-line, one-time only, data-driven
exploration phase. From this knowledge, Vito learns the optimal strategy
using RL (Sutton and Barto, 1998). The goal of Vito during the on-line
personalization phase is then to sequentially choose actions that maximize
future rewards, and therefore bring Vito to the state representing the solution of the personalization problem. To setup the algorithm, the user needs
to define what observations need to be matched, the allowed actions, and a
single hyper-parameter related to the desired granularity of the state-space.
Then everything is learned automatically. The algorithm does not depend
on the underlying model.
Vito was evaluated on three different tasks. First, in a synthetic experiment, convergence properties of the algorithm were analyzed. Then,
two tasks involving real clinical data were evaluated: the inverse problem
of cardiac electrophysiology and the personalization of a lumped-parameter
model of whole-body circulation. The obtained results suggested that Vito
can achieve equivalent (or better) goodness of fit as standard optimization
methods, increased robustness and faster convergence rates.
A number of novelties and improvements over Neumann et al. (2015)
are featured in this manuscript. First, an automatic, data-driven statespace quantization method is introduced that replaces the previous manual
technique. Second, the need to provide user-defined initial parameter values
is eliminated by employing a new data-driven technique to initialize personalization of unseen data. Third, a stochastic personalization policy is
introduced, for which the previously used standard deterministic policy is
a special case. Fourth, the convergence properties are evaluated in parameter space using a synthetic personalization scenario. In addition, thorough
evaluation of Vito’s performance with increasing amount of training samples was conducted and personalization of the whole-body circulation model
was extended to several variants involving two to six parameters. Finally,
the patient database used for experimentation was extended from 28 to 83
patients for the cardiac electrophysiology experiments, and from 27 to 56
for the whole-body circulation experiments.

4

The remainder of this manuscript is organized as follows. Sec. 2 presents
the method. In Sec. 3, the experiments are described and the results are
presented. Sec. 4 concludes the manuscript with a summary and discussions
about potential limitations and extensions of the method.
2. Method
This section presents the reinforcement-learning (RL) framework for
computational model personalization. Sec. 2.1 introduces Markov decision
process (MDP). Sec. 2.2 defines the personalization problem and how it can
be reformulated in terms of an MDP. Sec. 2.3 describes how the artificial
agent, Vito, learns how the model behaves. Next, Sec. 2.4 provides details
about state-space quantization, and Sec. 2.5 describes how the model knowledge is encoded in the form of transition probabilities. All steps mentioned
so far are performed in an off-line training phase. Finally, Sec. 2.6 explains
how the learned knowledge is applied on-line to personalize unseen data.
2.1. Model-based Reinforcement Learning
2.1.1. MDP Definition
A crucial prerequisite for applying RL is that the problem of interest, here personalization, can be modeled as a Markov decision process
(MDP). An MDP is a mathematical framework for modeling decision making when the decision outcome is partly random and partly controlled by
a decision maker (Sutton and Barto, 1998). Formally, an MDP is a tuple
M = (S, A, T , R, γ), where:
• S is the finite set of states that describe the agent’s environment, nS
is the number of states, and st ∈ S is the state at time t.
• A is the finite set of actions, which allow the agent to interact with
the environment, nA is the number of actions, and at ∈ A denotes the
action performed at time t.
• T : S × A × S → [0; 1] is the stochastic transition function, where
T (st , at , st+1 ) describes the probability of arriving in state st+1 after
the agent performed action at in state st .
• R : S × A × S → R is the scalar reward function, where rt+1 =
R(st , at , st+1 ) is the immediate reward the agent receives at time t + 1
after performing action at in state st resulting in state st+1 .
• γ ∈ [0; 1] is the discount factor that controls the importance of future
versus immediate rewards.

5

2.1.2. Value Iteration
The value of a state, V ∗ (s), is the expected discounted reward the agent
accumulates when it starts in state s and acts optimally in each step:
(∞
)
X
∗
k
V (s) = E
γ rt+k+1 st = s ,
(1)
k=0

where E{} denotes the expected value given the agent always selects the
optimal action, and t is any time step. Note that the discount factor γ is
a constant and the superscript k its exponent. V ∗ can be computed using
value iteration (Sutton and Barto, 1998), an iterative algorithm based on
dynamic programming. In the first iteration i = 0, let Vi : S → R denote
an initial guess for the value function that maps states to arbitrary values.
Further, let Qi : S × A → R denote the ith “state-action value function”guess, which is computed as:
X


Qi (s, a) =
T (s, a, s0 ) R(s, a, s0 ) + γVi (s0 ) .
(2)
s0 ∈S

Value iteration iteratively updates Vi+1 from the previous Qi :
∀s ∈ S :

Vi+1 (s) = max Qi (s, a) ,
a∈A

(3)

until the left- and right-hand side of Eq. 3 are equal for all s ∈ S; then
V ∗ ← Vi+1 and Q∗ ← Qi+1 . From this equality relation, also known as
the Bellman equation (Bellman, 1957), one can obtain an optimal problemsolving strategy for the problem described by the MDP (assuming that all
components of the MDP are known precisely). It is encoded in terms of a
deterministic optimal policy π ∗ : S → A:
π ∗ (s) = arg max Q∗ (s, a) ,

(4)

a∈A

i.e. a mapping that tells the agent in each state the optimal action to take.
2.1.3. Stochastic Policy
In this work not all components of the MDP are known precisely, instead some are approximated from training data. Value iteration, however,
assumes an exact MDP to guarantee optimality of the computed policy.
Therefore, instead of relying on the deterministic policy π ∗ (Eq. 4), a generalization to stochastic policies π̃ ∗ is proposed here to mitigate potential
issues due to approximations. Contrary to Eq. 4, where for each state only
the one action with maximum Q∗ -value is considered, a stochastic policy
stores several candidate actions with similar high Q∗ -value and returns one

6

Model
parameters

Model

Model
state

Objectives

Measured
data

Patient

Figure 2: A computational model f is a dynamic system that maps model input parameters x to model state (output) variables y. The goal of personalization is to tune x such
that the objectives c, defined as the misfit between y and the corresponding measured
data z of a given patient, are optimized (the misfit is minimized).

of them through a random process each time it is queried. To this end, the
Q∗ (s, ·)-values for a given state s are first normalized:
Q̃∗s (a) =

Q∗ (s, a) − mina0 ∈A [Q∗ (s, a0 )]
.
maxa0 ∈A [Q∗ (s, a0 )] − mina0 ∈A [Q∗ (s, a0 )]

(5)

All actions whose normalized Q̃∗s -value is below a threshold of  = 54 (set
empirically and used throughout the entire manuscript) are discarded, while
actions with large values are stored as potential candidates. Each time
the stochastic policy is queried, a = π̃∗ (s), it returns one of the candidate
actions P
a selected randomly with probability proportional to its Q̃∗s -value:
∗
Q̃s (a)/ a0 Q̃∗s (a0 ); the sum is over all candidate actions a0 .
2.2. Reformulation of the Model Personalization Problem into an MDP
2.2.1. Problem Definition
As illustrated in Fig. 2, any computational model f is governed by a set of
parameters x = (x1 , . . . , xnx )> , where nx denotes the number of parameters.
x is bounded within a physiologically plausible domain Ω, and characterized
by a number of ny (observable) state variables y = (y1 , . . . , yny )> . The
state variables can be used to estimate x. Note that some parameters may
be pre-estimated or assigned fixed values. The goal of personalization is to
optimize a set of nc objectives c = (c1 , . . . , cnc )> . The objectives are scalars
defined as ci = d(yi , zi ), where d is a measure of misfit, and zi denotes
the patient’s measured data (z) corresponding to yi . In this work d(yi , zi ) =
yi −zi . Personalization is considered successful if all user-defined convergence
criteria ψ = (ψ1 , . . . , ψnc )> are met. The criteria are defined in terms of
maximum acceptable misfit per objective: ∀i ∈ {1, . . . , nc } : |ci | < ψi .
2.2.2. Problem Reformulation
Personalization is mapped to a Markov decision process as follows:
States: An MDP state encodes the misfit between the computed model
state (outcome of forward model run) and the patient’s measurements.
Thus, MDP states carry the same type of information as objective vectors c,
yet the number of MDP states has to be finite (Sec. 2.1), while there are an
infinite number of different objective vectors due to their continuous nature.
Therefore the space of objective vectors in Rnc is reduced to a finite set of
7

representative states: the MDP states S, each s ∈ S covering a small region
in that space. One of those states, ŝ ∈ S, encodes personalization success
as it is designed such that it covers exactly the region where all convergence
criteria are satisfied. The goal of Vito is to learn how to reach that state.
Actions: Vito’s actions modify the parameters x to fulfill the objectives c.
An action a ∈ A consists in either in- or decrementing one parameter xi by
1×, 10× or 100× a user-specified reference value δi with δ = (δ1 , . . . , δnx )> .
This empirically defined quantization of the intrinsically continuous action
space yielded good results for the problems considered in this work.
Transition function: T encodes the agents knowledge about the computational model f and is learned automatically as described in Sec. 2.5.
Rewards: Inspired by the “mountain car” benchmark (Sutton and Barto,
1998), the rewards are defined as always being equal to R(s, a, s0 ) = −1
(punishment), except when the agent performs an action resulting in personalization success, i.e. when s0 = ŝ. In that case, R(·, ·, ŝ) = 0 (no punishment).
Discount factor: The large discount factor γ = 0.99 encourages policies
that favor future over immediate rewards, as Vito should always prefer the
long-term goal of successful personalization to short-term appealing actions
in order to reduce the risk of local minima.
2.3. Learning Model Behavior through Model Exploration
Like a human operator, Vito first learns how the model “behaves” by
experimenting with it. This is done through a “self-guided sensitivity analysis”. A batch of sample transitions is collected through model exploration
episodes E p = {ep1 , ep2 , . . . }. An episode epi is a sequence of ne-steps consecutive transitions generated from the model f and the patient p for whom
the target measurements zp are known. An episode is initiated at time
t = 0 by generating random initial model parameters xt within the physiologically plausible domain Ω. From the outputs of a forward model run
yt = f (xt ), the misfits to the patient’s corresponding measurements are
computed, yielding the objectives vector ct = d(yt , zp ). Next, a random exploration policy πrand that selects an action according to a discrete uniform
probability distribution over the set of actions is employed. The obtained
at ∈ A is then applied to the current parameter vector, yielding modified
parameter values xt+1 = at (xt ). From the output of the forward model run
yt+1 = f (xt+1 ) the next objectives ct+1 are computed. The next action at+1
is then selected according to πrand , and this process is repeated ne-steps − 1
times. Hence, each episode can be seen as a set of consecutive tuples:
e = {(xt , yt , ct , at , xt+1 , yt+1 , ct+1 ),

t = 0, . . . , ne-steps − 1} .

(6)

In this work, ne-steps = 100 transitions are created in each episode as a tradeoff between sufficient length of an episode to cover a real personalization
scenario and sufficient exploration of the parameter space.
8

The model is explored with many different training patients and the
S resulting episodes are combined into one large training episode set E = p E p .
The underlying hypothesis (verified in experiments) is that the combined
E allows to cancel out peculiarities of individual patients, i.e. to abstract
from patient-specific to model-specific knowledge.
2.4. From Computed Objectives to Representative MDP State
As mentioned above, the continuous space of objective vectors is quantized into a finite set of representative MDP states S. A data-driven approach is proposed. First, all objective vectors observed during training are
clustered according to their distance to each other. Because the ranges of
possible values for the individual objectives can vary significantly depending on the selected measurements (due to different types of measurements,
different units, etc.), the objectives should be normalized during clustering to avoid bias towards objectives with relatively large typical values. In
this work the distance measure performs implicit normalization to account
for these differences: the distance between two objective vectors (c1 , c2 ) is
defined relative to the inverse of the convergence criteria ψ:
q
kc1 − c2 kψ = (c1 − c2 )> diag(ψ)−1 (c1 − c2 ) ,
(7)
where diag(ψ)−1 denotes a diagonal matrix with ( ψ11 , ψ12 , . . . ) along its diagonal. The centroid of a cluster is the centroid of a representative state.
In addition, a special “success state” ŝ representing personalization success
is created, which covers the region in state-space where all objectives are
met: ∀i : |ci | < ψi . The full algorithm is described in Appendix A. Finally, an operator φ : Rnc → S that maps continuous objective vectors c to
representative MDP states is introduced:
φ(c) = arg min kc − ξs kψ

(8)

s∈S

where ξs denotes the centroid corresponding to state s. For an example
state-space quantization see Fig. 3.
2.5. Transition Function as Probabilistic Model Representation
In this work, the stochastic MDP transition function T encodes the
agent’s knowledge about the computational model f . It is learned from the
training data E. First, the individual samples (xt , yt , ct , at , xt+1 , yt+1 , ct+1 )
are converted to state-action-state transition tuples Ê = {(s, a, s0 )}, where
s = φ(ct ), a = at and s0 = φ(ct+1 ). Then, T is approximated from statistics
over the observed transition samples:
T (s, a, s0 ) = P

{(s, a, s0 ) ∈ Ê}
,

s00 ∈S

9

{(s, a, s00 ) ∈ Ê}

(9)

Data-driven quantization

Manual quantization

150

150

100

100

50

50

0

0

-50

-50

-100

-100

-150

-150

-150

-100

-50

0

50

100

150

-150

-100

-50

0

50

100

150

Figure 3: State-space quantization. Left: Example data-driven quantization of a twodimensional state-space into nS = 120 representative states. The states are distributed
according to the observed objective vectors c in one of the experiments in Sec. 3.2. The
objectives were QRS duration [ms] (c1 ) and electrical axis [deg] (c2 ). The center rectangle
(green region) denotes the success state ŝ where all objectives are met (∀i : |ci | < ψi ); see
text for details. Right: Manual quantization as used in Neumann et al. (2015).

where |{·}| denotes the cardinality of the set {·}. If nS and nA are large
compared to the total number of samples it may occur that some stateaction combinations are not observed: |{(s, a, ·) ∈ Ê}| = 0. In that case
uniformity is assumed: ∀s00 ∈ S : T (s, a, s00 ) = 1/nS .
M is now fully defined. Value iteration (Sec. 2.1) is applied and the
stochastic policy π̃∗ is computed, which completes the off-line phase.
2.6. On-line Model Personalization
On-line personalization, as illustrated in Fig. 4, can be seen as a two-step
procedure. First, Vito initializes the personalization of unseen patients from
training data. Second, Vito relies on the computed policy π̃∗ to guide the
personalization process.
2.6.1. Data-driven Initialization
Good initialization can be decisive for a successful personalization. Vito’s
strategy is to search for forward model runs in the training database E for
which the model state f (x) = y ≈ zp is similar to the patient’s measurements. To this end, Vito examines all parameters Ξ = {x ∈ E | f (x) ≈ zp }
that yielded model states similar to the patient’s measurements. Due to
ambiguities induced by the different training patients, data noise and model
assumptions, Ξ could contain significantly dissimilar parameters. Hence,
picking a single x ∈ Ξ might not yield the best initialization. Analyzing Ξ
probabilistically instead helps Vito to find likely initialization candidates.
The details of the initialization procedure are described in Appendix B.

10

Unseen
patient

Personalized
parameters

Initialization

Run forward model

Update parameters

Check convergence

Observe state

Select action

and detect oscillations

Figure 4: Vito’s probabilistic on-line personalization phase. See text for details.

Given the patient’s measurements zp , the procedure outputs a list of initialization candidates X0 = (x00 , x000 , . . . ). The list is sorted by likelihood with
the first element, x00 , being the most likely one.
2.6.2. Probabilistic Personalization
The first personalization step initializes the model parameter vector x0
with the most likely among all initialization candidates, x0 ∈ X0 (see previous section for details). Then, as illustrated in Fig. 4, Vito computes the
forward model y0 = f (x0 ) and the misfit between the model output and the
patient’s measurements c0 = d(y0 , zp ) to derive the first state s0 = φ(c0 ).
Given s0 , Vito decides from its policy the first action to take a0 = π̃∗ (s0 ),
and walks through state-action-state sequences to personalize the computational model f by iteratively updating the model parameters through MDP
actions. Bad initialization could lead to oscillations between states as observed in previous RL works (Kveton and Theocharous, 2012; Neumann
et al., 2015). Therefore, upon detection of an oscillation, which is done by
monitoring the parameter traces to detect recurring sets of parameter values,
the personalization is re-initialized at the second-most-likely x0 ∈ X0 , etc. If
all |X0 | initialization candidates have been tested, a potential re-initialization
defaults to fully random within the physiologically plausible parameter domain Ω. The process terminates once Vito reaches state ŝ (success), or when
a pre-defined maximum number of iterations is reached (failure).
3. Experiments
Vito was applied to a synthetic parameter estimation problem and to two
challenging problems involving real clinical data: personalization of a cardiac electrophysiology (EP), and a whole-body-circulation (WBC) model.
All experiments were conducted using leave-one-out cross-validation. The
numbers of datasets and transition samples used for the different experiments are denoted ndatasets and nsamples , respectively.
3.1. Synthetic Experiment: the Rosenbrock Function
First, Vito was employed in a synthetic scenario, where the ground-truth
model parameters were known. The goals were to test the ability of Vito to
11

5

≥5
0.25

4.5
4

2

3.5

0.25

3
1

3

0

2.5
2

0.25

-1
-2

1.5
1

-4

0.5

-5
-5 -4 -3 -2 -1

0.25

-3

0

1

2

3

4

5 -5 -4 -3 -2 -1

0

1

2

3

4

5 -5 -4 -3 -2 -1

0

1

2

3

4

5

Maximum observerd error

4

0

Figure 5: Synthetic experiment. Left: Contour plot of the Rosenbrock function f α=1 with
global minimum at x = (1, 1)> (red dot). The color scale is logarithmic for visualization
purposes: the darker, the lower the function value. Mid: Maximum L2 -error in parameter
space after personalization over all functions for varying initial parameter values. See text
for details. Yellow represents errors ≥ 5 (maximum observed error ≈ 110). Right: Same
as mid panel, except the extended action set was used. The red dots are the 100 groundtruth parameters x = (α, α2 )> generated for random α.

optimize cost functions generically, and to directly evaluate the performance
in the parameter space.
3.1.1. Forward Model Description
The Rosenbrock function (Rosenbrock, 1960), see Fig. 5, left panel, is
a non-convex function that is often used to benchmark optimization algorithms. It was treated as the forward model in this experiment:
f α (x1 , x2 ) = (α − x1 )2 + 100 · (x2 − x21 )2 ,

(10)

where x = (x1 , x2 )> were the model parameters to estimate for any α, and
f α : Ω → R. As described in Sec. 2.2.2, each of Vito’s actions a ∈ A in- or
decrements a parameter value by multiples (1×, 10×, 100×) of parameterspecific reference values. The reference values were set to δ = (0.01, 0.01)> ,
determined as 0.1% of the defined admissible parameter space per dimension,
Ω = [−5; 5]2 . The parameter α ∈ R defines a family of functions {f α }. The
goal was to find generically arg minx1 ,x2 f α (x1 , x2 ).
The Rosenbrock function has a unique global minimum at x = (α, α2 )> ,
where both terms T1 = (α − x1 ) and T2 = (x2 − x21 ) evaluate to 0. The
personalization objectives were therefore defined as c = (|T1 − 0|, |T2 − 0|)> ,
with the measured data z = (0, 0)> were zero for both objectives and the
computed data y = (T1 , T2 )> . The convergence criteria were set empirically
to ψ = (0.05, 0.05)> .
3.1.2. Evaluation
Vito was evaluated on ndatasets = 100 functions f α with randomly generated α ∈ [−2, 2]. In the off-line phase, for each function, nsamples =
12

10 · ne-steps = 1000 samples, i.e. ten training episodes, each consisting in
ne-steps = 100 transitions (Sec. 2.3), were generated to learn the policy. The
number of representative states was set to nS = 100. To focus on Vito’s
on-line personalization capabilities, both the data-driven initialization and
the re-initialization on oscillation (Sec. 2.6) were disabled. In total, 441 experiments with different initializations (sampled on a 21 × 21 uniform grid
spanned in Ω) were conducted. For each experiment all 100 functions were
personalized using leave-one-family-function-out cross validation, and the
error value from the function exhibiting the maximum L2 -error (worst-case
scenario) between ground-truth (α, α2 ) and estimated parameters was plotted. As one can see from the large blue region in Fig. 5, mid panel, for the
majority of initial parameter values Vito always converged to the solution
(maximum L2 -error < 0.25; the maximum achievable accuracy depended on
the specified convergence criteria ψ and on the reference values δ, which
“discretized” the parameter space). However, especially for initializations
far from the ground-truth (near border regions of Ω), Vito was unable to personalize some functions properly, which was likely due to the high similarity
of the Rosenbrock function shape in these regions.
To investigate this issue, the experiment was repeated after additional
larger parameter steps were added to the set of available actions: A0 =
A ∪ {±500δ1 ; ±500δ2 }. As shown in Fig. 5, right panel, Vito could now
personalize successfully starting from any point in Ω. The single spot with
larger maximum error (bright spot at approximately x = (−1, 2)> ) can be
explained by Vito’s stochastic behavior: Vito may have become unlucky if
it selected many unfavorable actions in sequence due to the randomness
introduced by the stochastic policy. Enabling re-initialization on oscillation
solved this issue entirely. In conclusion, this experiment showed that Vito
can learn how to minimize a cost function generically.
3.2. Personalization of Cardiac Electrophysiology Model
Vito was then tested in a scenario involving a complex model of cardiac
electrophysiology coupled with 12-lead ECG. Personalization was performed
for real patients from actual clinical data. A total of ndatasets = 83 patients were available for experimentation. For each patient, the end-diastolic
bi-ventricular anatomy was segmented from short-axis cine magnetic resonance imaging (MRI) stacks as described in Zheng et al. (2008). A tetrahedral anatomical model including myofibers was estimated and a torso atlas
affinely registered to the patient based on MRI scout images. See Zettinig
et al. (2014) for more details.
3.2.1. Forward Model Description
The depolarization time at each node of the tetrahedral anatomical
model was computed using a shortest-path graph-based algorithm, similar

13

to the one proposed in Wallman et al. (2012). Tissue anisotropy was modeled by modifying the edge costs to take into account fiber orientation. A
time-varying voltage map was then derived according to the depolarization
time: at a given time t, mesh nodes whose depolarization time was higher
than t were assigned a trans-membrane potential of −70 mV, 30 mV otherwise. The time-varying potentials were then propagated to a torso model
where 12-lead ECG acquisition was simulated, and QRS duration (QRSd)
and electrical axis (EA) were derived (Zettinig et al., 2014). The model
was controlled by the conduction velocities (in m/s) of myocardial tissue
and left and right Purkinje network: x = (vMyo , vLV , vRV )> . The latter two
domains were modeled as fast endocardial conducting tissue. The admissible parameter space Ω was set to [200; 1000] for vMyo and [500; 5000] for
both vLV and vRV . Reference increment values to build the action set A
were set to δ = (5, 5, 5)> m/s for the three model parameters. The goal of
EP personalization was to estimate x from the measured QRSd and EA.
Accounting for uncertainty in the measurements and errors in the model,
a patient was considered personalized if QRSd and EA misfits were below
ψ = (5 ms, 10◦ )> , respectively.
3.2.2. Number of Representative States
In contrast to Neumann et al. (2015), where state-space quantization
required manual tuning of various threshold values, the proposed approach
relies on a single hyper-parameter only: nS , the number of representative
states (Sec. 2.4). To specify nS , eight patients were selected for scouting.
Exhaustive search was performed for nS ∈ {10, 20, . . . , 490, 500} representative states. The goodness of a given configuration was evaluated based
on the success rate (relative number of successfully personalized cases according to convergence criteria ψ) over five independent, consecutive, leaveone-patient-out cross-validated personalization runs of the eight patients.
Furthermore, the average number of required forward model runs was considered. To this end, 100 training episodes (100·ne-steps = 104 transition samples) per patient were generated for each personalization run as described in
Sec. 2.3. As one can see from Fig. 6, good performance was achieved from 50
to 300 representative states. The large range of well performing nS indicates
a certain level of robustness with respect to that hyper-parameter. A slight
performance peak at 120 representative states was observed. Therefore,
nS = 120 was selected for further experimentation as compromise between
maintaining a low number of states and sufficient state granularity. An example quantization with nS = 120 is visualized in Fig. 3. The eight scouting
datasets were discarded for the following experiments to avoid bias in the
analysis.

14

0.8

80

0.6

60

0.4

40

0.2

20

0

10

50

100

150

200
250
300
350
Number of representative states

400

450

0
500

Average # forward model runs

100

120

Success rate

1

Figure 6: Hyper-parameter scouting. Vito’s performance for varying number of representative states nS on eight scouting datasets. The solid and dashed curves represent success
rate and average number of forward runs until convergence, respectively, aggregated over
five personalization runs with varying training data.

3.2.3. Reference Methods
Vito’s results were compared to two standard personalization methods
based on BOBYQA (Powell, 2009), a widely-used gradient-free optimizer
known for its robust performance and fast convergence. The first approach,
“BOBYQA simple”, mimicked the most basic estimation setup, where only
the minimum level of model and problem knowledge were assumed.
Pnc The
objective function was the sum of absolute QRSd and EA errors: i=1 |ci |.
It was minimized in a single optimizer run where all three parameters in
x were tuned simultaneously. The algorithm terminated once all convergence criteria ψ were satisfied (success) or if the number of forward model
evaluations exceeded 100 (failure). The second approach, “BOBYQA cascade”, implemented an advanced estimator with strong focus on robustness,
which computed the optimum parameters in a multi-step iterative fashion.
It is based on Seegerer et al. (2015), where tedious manual algorithm and
cost function tuning was performed on a subset of the data used in this
manuscript. In a first step, the myocardial conduction velocity was tuned to
yield good match between computed and measured QRS duration. Second,
left and right endocardial conduction velocities were optimized to minimize
electrical axis error. Both steps were repeated until all parameter estimates
were stable.
To remove bias towards the choice of initial parameter values, for each
of the two methods all datasets were personalized 100 times with different
random initializations within the range of physiologically plausible values
Ω. The differences in performance were striking: only by changing the
initialization, the number of successfully personalized cases varied from 13 to
37 for BOBYQA simple, and from 31 to 51 for BOBYQA cascade (variability
of more than 25% of the total number of patients). These results highlight
the non-convexity of the cost function to minimize.

15

0

135

Data-driven
initialization

90

Full personalization

Absolute EA error [deg]

50

Full personalization

Fixed initialization

Absolute QRSd error [ms]

100

Data-driven
initialization

Fixed initialization

180

150

45

0

Training samples per dataset

Training samples per dataset

Figure 7: Absolute errors over all patients after initialization with fixed parameter values
(blue), after data-driven initialization for increasing amount of training data (white), and
after full personalization with Vito (green). Data-driven initialization yielded significantly
reduced errors if sufficient training data were available (> 102 ) compared to initialization
with fixed values. Full personalization further reduced the errors by a significant margin.
The red bar and the box edges indicate the median absolute error, and the 25 and 75
percentiles, respectively. Left: QRS duration errors. Right: Electrical axis errors.

3.2.4. Full Personalization Performance
First, Vito’s overall performance was evaluated. The full personalization
pipeline consisting in off-line learning, initialization, and on-line personalization was run on all patients with leave-one-patient-out cross-validation
using 1000 training episodes (nsamples = 1000 · ne-steps = 105 transition samples) per patient. The maximum number of iterations was set to 100. The
green box plots in the two panels of Fig. 7 summarize the results. The
mean absolute errors were 4.1 ± 5.6 ms and 12.4 ± 13.3◦ in terms of QRSd
and EA, respectively, a significant improvement over the residual error after
initialization. In comparison to the reference methods, the best BOBYQA
simple run yielded absolute errors of 4.4 ± 10.8 ms QRSd and 15.5 ± 18.6◦
EA on average, and the best BOBYQA cascade run 0.1 ± 0.2 ms QRSd and
11.2 ± 15.8◦ EA, respectively. Thus, in terms of EA error all three methods
yielded comparable performance, and while BOBYQA simple and Vito performed similarly in terms of QRSd, BOBYQA cascade outperformed both
in this regard. However, considering success rates, i.e. successfully personalized patients according to the defined convergence criteria (ψ) divided by
total number of patients, both the performance of Vito (67%) and BOBYQA
cascade (68%) were equivalent, while BOBYQA simple reached only 49% or
less. In terms of run-time, i.e. average number of forward model runs until convergence, Vito (31.8) almost reached the high efficiency of BOBYQA
simple (best: 20.1 iterations) and clearly outperformed BOBYQA cascade
(best: 86.6 iterations), which means Vito was ≈ 2.5× faster.

16

3.2.5. Residual Error after Initialization
A major advantage over standard methods such as the two BOBYQA approaches is Vito’s automated, data-driven initialization method (Sec. 2.6.1),
which eliminates the need for user-provided initial parameter values. To
evaluate the utility of this step, personalization using Vito was stopped directly after initialization (the most likely x0 was used) and the errors in terms
of QRSd and EA resulting from a forward model run f with the computed
initial parameter values were quantified. This experiment was repeated for
increasing number of transition samples per dataset: nsamples = 100 . . . 105 ,
and the results were compared to the error after initialization when fixed
initial values were used (the initialization of the best performing BOBYQA
experiment was used). As one can see from Fig. 7, with increasing amount
of training data both errors decreased notably. As few as 102 transitions per
dataset already provided more accurate initialization than the best tested
fixed initial values. Thus, not only does this procedure simplify the setup
of Vito for new problems (no user-defined initialization needed), this experiment showed that it can reduce initial errors by a large margin, even
when only few training transitions were available. It should be noted that
Vito further improves the model fit in its normal operating mode (continue
personalization after initialization), as shown in the previous experiment.
3.2.6. Convergence Analysis
An important question in any RL application relates to the amount of
training needed until convergence of the artificial agent’s behavior. For Vito
in particular, this translates to the amount of transition samples required
to accurately estimate the MDP transition function T to compute a solid
policy on the one hand, and to have enough training data for reliable parameter initialization on the other hand. To this end, Vito’s overall performance
(off-line learning, initialization, personalization) was evaluated for varying
number of training transition samples per dataset. As one can see from the
results in Fig. 8, with increasing amount of training data the performance
increased, suggesting that the learning process was working properly. Even
with relatively limited training data of only nsamples = 102 samples per patient, Vito outperformed the best version of BOBYQA simple (49% success
rate). Starting from nsamples ≈ 3000, a plateau at ≈66% success rate was
reached, which remained approximately constant until the maximum tested
number of samples. This was almost on par with the top BOBYQA cascade
performance (68% success rate). Also the run-time performance increased
with more training data. For instance, Vito’s average number of iterations
was 36.2 at 103 samples, 31.5 at 104 samples, or 31.8 at 105 samples.
These results suggested that not only Vito can achieve similar performance as an advanced, manually engineered method, but also the number
of required training samples was not excessive. In fact, a rather limited
and thus well manageable amount of data, which can be computed in a
17

60

0.4

40
# Forward model runs

0.2
0
100

101

102
103
Training samples per dataset

104

20
0

0.6
0.4
0.2
0

BOBYQA
simple

BOBYQA
cascade

100
80
60
40
20
0

Average # forward model runs

0.6

# Forward model runs

0.8

Success rate

80

Reference Methods
# Forward model runs

Success rate

1

Success rate

Success rate

0.8

100

Average # forward model runs
Success rate

Vito

1

Figure 8: EP personalization results. Personalization success rate in blue and average
number of iterations in red. Left: Vito’s performance for increasing number of training
transition samples per dataset. Each dot represents results from one experiment (crossvalidated personalization of all 75 datasets), solid/dashed line is low-pass filtered mean,
shaded areas represent 0.5× and 1× standard deviation. Right: Performance of both
reference methods. Each shade represents 10% of the results, sorted by performance.

reasonable time-frame, sufficed.
3.3. Personalization of Whole-Body Circulation Model
Next, Vito was asked to personalize a lumped-parameter whole-body
circulation (WBC) model from pressure catheterization and volume data.
A subset of ndatasets = 56 patients from the EP experiments were used for
experimentation. The discrepancy was due to missing catheterization data
for some patients, which was required for WBC personalization only. For
each patient, the bi-ventricular anatomy was segmented and tracked from
short-axis cine MRI stacks throughout one full heart cycle using shapeconstraints, learned motion models and diffeomorphic registration (Wang
et al., 2013). From the time-varying endocardial meshes, ventricular volume
curves were derived. Manual editing was performed whenever necessary.
3.3.1. Forward Model Description
The WBC model to personalize was based on Itu et al. (2014). It contained a heart model (left ventricle (LV) and atrium, right ventricle and
atrium, valves), the systemic circulation (arteries, capillaries, veins) and the
pulmonary circulation (arteries, capillaries, veins). Time-varying elastance
models were used for all four chambers of the heart. The valves were modeled
through a resistance and an inertance. A three-element Windkessel model
was used for the systemic and pulmonary arterial circulation, while a twoelement Windkessel model was used for the systemic and pulmonary venous
circulation. We refer the reader to Itu et al. (2014); Neumann et al. (2015);
Westerhof et al. (1971) for more details. Personalization was performed with
respect to the patient’s heart rate as measured during catheterization.
The goal of this experiment was to compare Vito’s personalization performance for the systemic part of the model in setups with increasing number of parameters to tune and objectives to match. To this end, Vito was
18

x
Initial volume
LV max. elastance
Aortic resistance
Aortic compliance
Dead volume
Time to Emax

Default value
400 mL
2.4 mmHg/mL
1100 g/(cm4 s)
1.4 ·109 cm4 s2 /g
10 mL
300 ms

Ω
[200; 1000] mL
[0.2; 5] mmHg/mL
[500; 2500] g/(cm4 s)
[0.5; 6] ·109 cm4 s2 /g
[−50; 500] mL
[100; 600] ms

Setups
6, 5, 3, 2
6, 5, 3, 2
6, 5, 3
6, 5
6, 5
6

Table 1: WBC parameters x, their default values and domain Ω. The last column denotes the experiment setups in which a parameter was personalized (e.g. “5”: parameter
was among the estimated parameters in 5p experiment). Default values were used in
experiments where the respective parameters were not personalized.

c
End-diastolic LV volume
End-systolic LV volume
Mean aortic pressure
Peak-systolic aortic pressure
End-diastolic aortic pressure
Ejection time

ψ
20 mL
20 mL
10 mmHg
10 mmHg
10 mmHg
50 ms

Measured range
[129; 647] mL
[63; 529] mL
[68; 121] mmHg
[83; 182] mmHg
[48; 99] mmHg
[115; 514] ms

Setups
6, 5, 3, 2
6, 5, 3, 2
6, 5, 3
6, 5
6, 5
6

Table 2: WBC objectives c, their convergence criteria ψ and range of measured values in
the patient population used for experimentation.

employed on setups with two to six parameters (2p, 3p, 5p, 6p): initial
blood volume, LV maximum elastance, time until maximum elastance is
reached, total aortic resistance and compliance, and LV dead volume. The
reference values δ to define Vito’s allowed actions A were set to .5% of the
admissible parameter range Ω for each individual parameter, see Table 1 for
details. The personalization objectives were MRI-derived end-diastolic and
end-systolic LV volume, ejection time (time duration during which the aortic valve is open and blood is ejected), and peak-systolic, end-diastolic, and
mean aortic blood pressures as measured during cardiac catheterization, see
Fig. 9. To account for measurement noise, personalization was considered
successful if the misfits per objective were below acceptable threshold values
ψ as listed in Table 2.
3.3.2. Number of Representative States
Along the same lines as Sec. 3.2.2, the hyper-parameter for state-space
quantization was tuned based on the eight scouting patients. The larger the
dimensionality of the state-space, the more representative states were needed
to yield good performance. In particular, for the different WBC setups, the
numbers of representative states (nS ) yielding the best scouting performance
were 70, 150, 400 and 600 for the 2p, 3p, 5p and 6p setup, respectively. The
scouting datasets were discarded for the following experiments.
19

LV volume [mL]

Aortic pressure [mmHg]

160
140

WBC 2p
Measured
Computed

160

WBC 3p

160

140

140

120

120

120

100

100

100

80

80

60

60

60

40

40

40

20
200

20
200

20
200

185

End-diastolic
volume

185

Mean
pressure

80

End-diastolic
volume

185

WBC 5p
Peak-systolic
pressure

160
140
100

Mean
pressure

80

End-diastolic
volume

40
20
200

End-diastolic
volume
Ejection
time

170

155

155

155

155

140

140

140

110

0 0.2 0.4 0.6 0.8 1

Time [s]

125
110

125
110

0 0.2 0.4 0.6 0.8 1

Time [s]

Endsystolic
volume
0 0.2 0.4 0.6 0.8 1

Time [s]

End-diastolic
aortic pressure

185

170

125

Mean
pressure

60
End-diastolic
aortic pressure

170

Endsystolic
volume

Peak-systolic
pressure

120

170

Endsystolic
volume

WBC 6p

140
125
110

Endsystolic
volume
0 0.2 0.4 0.6 0.8 1

Time [s]

Figure 9: Goodness of fit in terms of time-varying LV volume and aortic pressure for Vito
personalizing an example patient based on the different WBC setups. The added objectives
per setup are highlighted in the respective column. With increasing number of parameters
and objectives Vito manages to improve the fit between model and measurements.

3.3.3. Reference Method
A gradient-free optimizer (Lagarias et al., 1998) based on the simplex
method was used to benchmark Vito. The objective function was the sum
of squared differences between computed and measured values, weighted by
the inverse of the convergence criteria to counter the different ranges of
objective values (e.g. due to different types of measurements and different
units): kckψ (Eq. 7). Compared to non-normalized optimization, the algorithm converged up to 20% faster and success rates increased by up to
8% under otherwise identical conditions. Personalization was terminated
once all convergence criteria were satisfied (success), or when the maximum
number of iterations was reached (failure). To account for the increasing
complexity of optimization with increasing number of parameters nx , the
maximum number of iterations was set to 50 · nx for the different setups.
As one can see from Fig. 10, right panels, with increasing number of
parameters to be estimated, the performance in terms of success rate and
number of forward model runs decreased slightly. This is expected as the
problem becomes harder. To suppress bias originating from (potentially
poor) initialization, the reference method was run 100 times per setup (as
in EP experiments), each time with a different, randomly generated set of
initial parameter values. The individual performances varied significantly
for all setups.

20

Vito

Success rate

1

1

0.8

0.8

0.6

0.6

0.4

0.4
2p
3p
5p
6p

Average # forward model runs

0.2
0
0
10
300

1

10

2

10

3

0.2
0

4

10

10

2p

3p

5p

6p

300
2p
3p
5p
6p

250
200

250
200

150

150

100

100

50

50

0
100

Reference Method

101

102
103
Training samples per dataset

104

0

2p

3p
5p
#Parameters

6p

Figure 10: WBC model personalization results (top: success rate, bottom: average number
of forward model runs until convergence) for various estimation setups (different colors),
see text for details. Left: Vito’s performance for increasing number of training transition
samples per dataset. Each dot represents results from one experiment (cross-validated
personalization of all 48 datasets), solid/dashed lines are low-pass filtered mean, shaded
areas represent 0.5× and 1× standard deviation. Right: Performance of reference method.
Each shade represents 10% of the results, sorted by performance; darkest shade: best 10%.

21

3.3.4. Convergence Analysis
For each WBC setup the full Vito personalization pipeline was evaluated
for increasing training data (nsamples = 100 . . . 105 ) using leave-one-patientout cross-validation. The same iteration limits as for the reference method
were used. The results are presented in Fig. 10, left panels. With increasing
data, Vito’s performance, both in terms of success rate and run-time (iterations until convergence), increased steadily until reaching a plateau. As
one would expect, the more complex the problem, i.e. the more parameters
and objectives involved in the personalization, the more training data was
needed to reach the same level of performance. For instance, Vito reached
80% success rate with less than nsamples = 50 training samples per dataset in
the 2p setup, whereas almost 90× as many samples were required to achieve
the same performance in the 6p setup.
Compared to the reference method, given enough training data, Vito
reached equivalent or better success rates (e.g. up to 11% higher success
rate for 6p) while significantly outperforming the reference method in terms
of run-time performance. In the most basic setup (2p), if nsamples ≥ 103 , Vito
converged after 3.0 iterations on average, while the best reference method
run required 22.6 iterations on average, i.e. Vito was seven times faster. For
the more complex setups (3p, 5p, 6p), the speed-up was not as drastic. Yet,
in all cases Vito outperformed even the best run of the reference method by
a factor of 1.8 or larger.
4. Conclusion
4.1. Summary and Discussion
In this manuscript, a novel personalization approach called Vito has been
presented. To our knowledge, it is the first time that biophysical model
personalization is addressed using artificial intelligence concepts. Inspired
by how humans approach the personalization problem, Vito first learns the
characteristics of the computational model under consideration using a datadriven approach. This knowledge is then utilized to learn how to personalize
the model using reinforcement learning. Vito is generic in the sense that it
requires only minimal and intuitive user input (parameter ranges, authorized
actions, number of representative states) to learn by itself how to personalize
a model.
Vito was applied to a synthetic scenario and to two challenging personalization tasks in cardiac computational modeling. The problem setups
and hyper-parameter configurations are listed in Table 3. In most setups
the majority of hyper-parameters were identical and only few (nS ) required
manual tuning, suggesting good generalization properties of Vito. Another
key result was that Vito was up to 11% more robust (higher success rates)
compared to standard personalization methods. Vito’s ability to generalize

22

Application
Rosenbrock
Rosenbrock ext.
EP
WBC 2p
WBC 3p
WBC 5p
WBC 6p

nx
2
2
3
2
3
5
6

nc
2
2
2
2
3
5
6

ndatasets
100
100
83 (75)
56 (48)
56 (48)
56 (48)
56 (48)

nS
100
100
120
70
150
400
600

nA /nx
6
8
6
6
6
6
6

nplateau
n/a
n/a
3 000
450
2 000
3 500
20 000

Table 3: Applications considered in this manuscript described in terms of the number
of parameters (nx ), objectives (nc ) and datasets (ndatasets ) used for experimentation (in
brackets: excluding scouting patients, if applicable); and Vito’s hyper-parameters: the
number of representative MDP states (nS ) and the number of actions per parameter
(nA /nx ). The last column (nplateau ) denotes the approximate number of samples needed
to reach the performance “plateau” (see convergence analyses in Sec. 3.2.6 and Sec. 3.3.4).

the knowledge obtained from a set of training patients to personalize unseen patients was shown as all experiments reported in this manuscript were
based on cross-validation. Furthermore, Vito’s robustness against training
patients for whom we could not find a solution was tested. In particular, for
about 20% of the patients, in none of the electrophysiology experiments in
Sec. 3.2 any personalization (neither Vito nor the reference methods) could
produce a result that satisfied all convergence criteria. Hence, for some
patients no solution may exist under the given electrophysiology model configuration1 . Still, all patients were used to train Vito, and surprisingly Vito
was able to achieve equivalent success rate as the manually engineered personalization approach for cardiac EP.
Generating training data could be considered Vito’s computational bottleneck. However, training is i) performed off-line and one-time only, and
ii) it is independent for each training episode and each patient. Therefore,
large computing clusters could be employed to perform rapid training by
parallelizing this phase. On-line personalization, on the contrary, is not parallelizable in its current form: the parameters for each forward model run
depend on the outcome of the previous iteration. Since the forward computations are the same for every “standard” personalization method (not
including surrogate-based approaches), the number of forward model runs
until convergence was used for benchmarking: Vito was up to seven times
faster compared to the reference methods. The on-line overhead introduced
by Vito (convert data into an MDP state, then query policy) is negligible.
As such, Vito could become a unified framework for personalization of
any computational physiological model, potentially eliminating the need for
1
Potential solution non-existence may be due to possibly invalid assumptions of the
employed EP model for patients with complex pathologies.

23

an expert operator with in-depth knowledge to design and engineer complex
optimization procedures.
4.2. Challenges and Outlook
Important challenges still remain, such as the incorporation of continuous
actions, the definition of states and their quantization. In this work we propose a data-driven state-space quantization strategy. Contrary to Neumann
et al. (2015), where a threshold-based state-quantization involving several
manually tuned threshold values (Fig. 3) was employed, the new method
is based on a single hyper-parameter only: the number of representative
states. Although it simplifies the setup of Vito, this quantization strategy
may still not be optimal, especially if only little training data is available.
Therefore, advanced approaches for continuous reinforcement learning with
value function approximation (Mnih et al., 2015; Sutton and Barto, 1998)
could be integrated to fully circumvent quantization issues.
At the same time, such methods could improve Vito’s scalability towards
high-dimensional estimation tasks. In this work we showed that Vito can be
applied to typical problems emerging in cardiac modeling, which could be
described as medium-scale problems with moderate number of parameters to
personalize and objectives to match. In unreported experiments involving
>10 parameters, however, Vito could no longer reach satisfactory performance, which is likely due to the steeply increasing number of transition
samples needed to sample the continuous state-space of increasing dimensionality sufficiently during training. The trends in Sec. 3.3 confirm the need
for more data. In the future, experience replay (Adam et al., 2012; Lin, 1993)
or similar techniques could be employed to increase training data efficiency.
Furthermore, massively parallel approaches (Nair et al., 2015) are starting
to emerge, opening up new avenues for large-scale reinforcement learning.
Although the employed reinforcement learning techniques guarantee convergence to an optimal policy, the computed personalization strategy may
not be optimal for the model under consideration as the environment is only
partially observable and the personalization problem ill-posed: there is no
guarantee for solution existence or uniqueness. Yet, we showed that Vito
can solve personalization more robustly and more effectively than standard
methods under the same conditions. However, a theoretical analysis in terms
of convergence guarantees and general stability of the method would be desirable, in particular with regards to the proposed re-initialization strategy.
As a first step towards this goal, in preliminary (unreported) experiments on
the EP and the WBC model we observed that the number of patients which
do not require re-initialization (due to oscillation) to converge to a successful
personalization consistently increased with increasing training data.
The data-driven initialization proposed in this work simplifies Vito’s
setup by eliminating the need for user-provided initialization. However,
currently there is no guarantee that the first initialization candidate is the
24

one that will yield the “best” personalization outcome. Therefore, one could
investigate the benefits of a fuzzy personalization scheme: many personalization processes could be run in parallel starting from the different initialization candidates. Parameter uncertainty quantification techniques (Neumann
et al., 2014a) could then be applied to compute a probability density function over the space of model parameters. Such approaches aim to gather
complete information about the solution-space, which can be used to study
solution uniqueness and other interesting properties.
An important characteristic of any personalization algorithm is its stability against small variations of the measured data. A preliminary experiment indicated good stability of Vito: the computed parameters from
several personalization runs, each involving small random perturbations of
the measurements, were consistent. Yet in a small group of patients some
parameter variability was observed, however, it was below the variability of
the reference method under the same conditions. To what extent certain
degrees of variability will impact other properties of the personalized model
such as its predictive power will be subject of future research. We will also
investigate strategies to improve Vito’s stability further. For instance, the
granularity of the state-space could provide some flexibility to tune the stability: less representative states means a larger region in state space per
state, thus small variations in the measured data might have less impact
on personalization outcome. However, this could in turn have undesirable
effects on other properties of Vito such as success rate or convergence speed
(see Sec. 3.2.2).
Beyond these challenges, Vito showed promising performance and versatility, making it a first step towards an automated, self-taught model personalization agent. The next step will be to investigate the predictive power
of the personalized models, for instance for predicting acute or long-term
response in cardiac resynchronization therapy (Kayvanpour et al., 2015; Sermesant et al., 2009).
Appendix A. Data-driven State-Space Quantization
This section describes the details of the proposed data-driven quantization approach to define the set of representative MDP states S (see Sec. 2.4).
It is based on clustering, in particular on the weighted k -means algorithm
described in Arthur and Vassilvitskii (2007). To this end, all objective vectors C = {c ∈ E} are extracted from the training data (Sec. 2.3). C ⊂ Rnc
represents all observed “continuous states”. The goal is to convert C into
the finite set of representative MDP states S while taking into account that
Vito relies on a special “success state” ŝ encoding personalization success.
The success state ŝ does not depend on the data, but on the maximum
acceptable misfit ψ. In particular, since personalization success implies that
all objectives are met, ŝ should approximate a hyperrectangle centered at
25

Figure A.11: Preprocessing of k -means input data to enforce the success state ŝ. Left:
Continuous state-space with observed objective vectors c (blue points). The points with
dashed outline will be canceled out. Right: Delineation of ŝ in green, enforced by inserted
vectors (green / red points) with large weights. See text for details.

0 and bounded at ±ψ, i.e. a small region in Rnc where ∀i : |ci | < ψi . To
enforce ŝ, the input to weighted k -means is preprocessed as follows.
First, the 0-vector is inserted into C, along with two vectors per dimension i, where all components are zero, except the ith component, which is set
to ±2ψi . These 2nc + 1 inserted vectors are later converted into centroids
of representative states to delineate the desired hyperrectangle for ŝ as illustrated in Fig. A.11. Furthermore, to avoid malformation of ŝ, no other
representative state should emerge within that region. Therefore, all vectors
c ∈ C, where ∀i : |ci | < 2ψi (except for the inserted vectors) are canceled
out by assigning zero weight, while the inserted vectors are assigned large
weights → ∞ and all remaining vectors weights of 1.
Next, k -means is initialized by placing a subset of the initial centroids at
the locations of the inserted states, and the remaining nS − 2nc − 1 centroids
at random vectors in C. Both the large weight and the custom initialization
enforce the algorithm to converge to a solution where one cluster centroid
is located at each inserted vector, while the other centroids are distributed
according to the training data. To ensure equal contribution of all objectives
(cancel out different units, etc.), similarity is defined relative to the inverse
of the user-defined convergence criteria (Eq. 7).
Finally, after k -means converged, the resulting centroids, denoted ξs , are
used to delineate the region in Rnc assigned to a representative state s.
Appendix B. Data-driven Initialization
This section describes the details of the proposed data-driven initialization approach to compute a list of candidate initialization parameter vectors
X0 = (x00 , x000 , . . . ) for a new patient p based on the patient’s measurements
zp and the training database E (see Sec. 2.6.1).
26

First, all model states are extracted from the training database: Υ =
{y ∈ E}. Next, Υ is fed to a clustering algorithm (e.g. k -means). As in
Appendix A, the distance measure is defined relative to the inverse of the
convergence criteria (Eq. 7). The output is a set of centroids (for simplicity,
in this work the number of centroids was set to nS ), and each vector is
assigned to one cluster based on its closest centroid. Let Υp ⊆ Υ denote the
members of the cluster whose centroid is closest to zp and Ξp = {x ∈ E |
f (x) ∈ Υp } the set of corresponding model parameters. For each cluster, an
approximation of the likelihood over the generating parameters is computed
in terms of a probability density function. In this work a Gaussian mixture
model is assumed:
GMMp (x) =

M
X

νm N (x; µm , Σm ) .

(B.1)

m=1

The parameter vectors in Ξp are treated as random samples drawn from
GMMp . Its properties, namely the number of mixture components M , their
weights νm , and their means µm and covariance matrices Σm , are estimated
from these samples using a multivariate kernel density estimator with automated kernel bandwidth estimation, see Kristan et al. (2011) for more
details. Finally, the M estimated means are selected as initialization candidates and stored in a list X0 = (µm0 , µm00 , . . . ). The elements of X0 are
sorted in descending order according to their corresponding νm -values to prioritize more likely initializations: µm0 is the mean with m0 = arg maxm νm .
References
References
Adam, S., Busoniu, L., Babuska, R., 2012. Experience replay for real-time
reinforcement learning control. IEEE Sys. Man. Cybern. 42 (2), 201–212.
Aguado-Sierra, J., Kerckhoffs, R. C. P., Lionetti, F., Hunt, D., Villongco,
C., Gonzales, M., Campbell, S. G., McCulloch, A. D., 2010. A computational framework for patient-specific multi-scale cardiac modeling. In:
Kerckhoffs, R. C. (Ed.), Patient-Specific Modeling of the Cardiovascular
System. Springer, pp. 203–223.
Aguado-Sierra, J., Krishnamurthy, A., Villongco, C., Chuang, J., Howard,
E., Gonzales, M. J., Omens, J., Krummen, D. E., Narayan, S., Kerckhoffs,
R. C. P., McCulloch, A. D., 2011. Patient-specific modeling of dyssynchronous heart failure: a case study. Prog. Biophys. Mol. Bio. 107 (1),
147–155.
Arthur, D., Vassilvitskii, S., 2007. k-means++: The advantages of careful
seeding. In: ACM-SIAM Symp. Discr. Algorithm. pp. 1027–1035.
27

Augenstein, K. F., Cowan, B. R., LeGrice, I. J., Nielsen, P. M., Young, A. A.,
2005. Method and apparatus for soft tissue material parameter estimation
using tissue tagged magnetic resonance imaging. J. Biomech. Eng. 127 (1),
148–157.
Barreto, A., Precup, D., Pineau, J., 2014. Practical kernel-based reinforcement learning. arXiv preprint arXiv:1407.5358.
Bellman, R., 1957. Dynamic Programming. Princeton University Press.
Bishop, C. M., 2006. Pattern recognition and machine learning. Vol. 4.
Springer New York.
Chabiniok, R., Moireau, P., Lesault, P.-F., Rahmouni, A., Deux, J.-F.,
Chapelle, D., 2012. Estimation of tissue contractility from cardiac cinemri using a biomechanical heart model. Biomech. Model. Mechan. 11 (5),
609–630.
Clayton, R., Bernus, O., Cherry, E., Dierckx, H., Fenton, F., Mirabella,
L., Panfilov, A., Sachse, F., Seemann, G., Zhang, H., 2011. Models of
cardiac tissue electrophysiology: progress, challenges and open questions.
Prog. Biophys. Mol. Bio. 104 (1), 22–48.
Delingette, H., Billet, F., Wong, K. C. L., Sermesant, M., Rhode, K., Ginks,
M., Rinaldi, C. A., Razavi, R., Ayache, N., 2012. Personalization of cardiac motion and contractility from images using variational data assimilation. IEEE T. Biomed. Eng. 59 (1), 20–24.
Frangi, A. F., Niessen, W. J., Viergever, M., 2001. Three-dimensional modeling for functional analysis of cardiac images, a review. IEEE T. Med. Imaging 20 (1), 2–5.
Hunter, P. J., Borg, T. K., 2003. Integration from proteins to organs: the
physiome project. Nat. Rev. Mol. Cell Bio. 4 (3), 237–243.
Itu, L., Sharma, P., Georgescu, B., Kamen, A., Suciu, C., Comaniciu, D.,
2014. Model based non-invasive estimation of PV loop from echocardiography. IEEE Eng. Med. Biol. Soc., 6774–6777.
Kaelbling, L. P., Littman, M. L., Moore, A. W., 1996. Reinforcement learning: A survey. J. Artif. Intell. Res., 237–285.
Kayvanpour, E., Mansi, T., Sedaghat-Hamedani, F., Amr, A., Neumann, D.,
Georgescu, B., Seegerer, P., Kamen, A., Haas, J., Frese, K. S., Irawati, M.,
Wirsz, E., King, V., Buss, S., Mereles, D., Zitron, E., Keller, A., Katus,
H. A., Comaniciu, D., Meder, B., 2015. Towards personalized cardiology:
Multi-scale modeling of the failing heart. PLoS ONE 10 (7), e0134869.

28

Kerckhoffs, R. C. P., Lumens, J., Vernooy, K., Omens, J., Mulligan, L., Delhaas, T., Arts, T., McCulloch, A., Prinzen, F., 2008. Cardiac resynchronization: insight from experimental and computational models. Prog. Biophys. Mol. Bio. 97 (2), 543–561.
Konukoglu, E., Relan, J., Cilingir, U., Menze, B. H., Chinchapatnam, P.,
Jadidi, A., Cochet, H., Hocini, M., Delingette, H., Jaı̈s, P., Haı̈ssaguerre,
M., Ayache, N., Sermesant, M., 2011. Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to eikonal-diffusion models in cardiac electrophysiology. Prog. Biophys. Mol. Bio. 107 (1), 134–146.
Krishnamurthy, A., Villongco, C. T., Chuang, J., Frank, L. R., Nigam, V.,
Belezzuoli, E., Stark, P., Krummen, D. E., Narayan, S., Omens, J. H.,
McCulloch, A. D., Kerckhoffs, R. C. P., 2013. Patient-specific models of
cardiac biomechanics. J. Comput. Phys. 244, 4–21.
Kristan, M., Leonardis, A., Skočaj, D., 2011. Multivariate online kernel
density estimation with gaussian kernels. Pattern Recogn. 44 (10), 2630–
2642.
Kuijpers, N. H., Hermeling, E., Bovendeerd, P. H., Delhaas, T., Prinzen,
F. W., 2012. Modeling cardiac electromechanics and mechanoelectrical
coupling in dyssynchronous and failing hearts. J. Cardiovasc. Transl. Res.
5 (2), 159–169.
Kveton, B., Theocharous, G., 2012. Kernel-based reinforcement learning on
representative states. In: Association for the Advancement of Artificial
Intelligence. pp. 977–983.
Lagarias, J. C., Reeds, J. A., Wright, M. H., Wright, P. E., 1998. Convergence properties of the Nelder-Mead simplex method in low dimensions.
SIAM J. Optimiz. 9 (1), 112–147.
Le Folgoc, L., Delingette, H., Criminisi, A., Ayache, N., 2013. Current-based
4D shape analysis for the mechanical personalization of heart models. In:
Medical Computer Vision. Recognition Techniques and Applications in
Medical Imaging. Vol. 7766 of LNCS. Springer, pp. 283–292.
Lin, L.-J., 1993. Reinforcement learning for robots using neural networks.
Tech. rep., DTIC Document.
Marchesseau, S., Delingette, H., Sermesant, M., Cabrera-Lozoya, R., TobonGomez, C., Moireau, P., Figueras i Ventura, R. M., Lekadir, K., Hernandez, A., Garreau, M., Donal, E., Leclercq, C., Duckett, S. G., Rhode,
K., Rinaldi, C. A., Frangi, A. F., Razavi, R., Chapelle, D., Ayache, N.,
2013. Personalization of a cardiac electromechanical model using reduced
29

order unscented kalman filtering from regional volumes. Med. Image Anal.
17 (7), 816–829.
Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare,
M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., Petersen, S., Beattie, C., Sadik, A., Antonoglou, I., King, H., Kumara, D.,
Wierstra, D., Legg, S., Hassabis, D., 2015. Human-level control through
deep reinforcement learning. Nature 518 (7540), 529–533.
Mülling, K., Kober, J., Kroemer, O., Peters, J., 2013. Learning to select and
generalize striking movements in robot table tennis. Int. J. Robot. Res.
32 (3), 263–279.
Nair, A., Srinivasan, P., Blackwell, S., Alcicek, C., Fearon, R., De Maria,
A., Panneershelvam, V., Suleyman, M., Beattie, C., Petersen, S., Legg, S.,
Mnih, V., Kavukcuoglu, K., Silver, D., 2015. Massively parallel methods
for deep reinforcement learning. arXiv:1507.04296.
Neumann, D., Mansi, T., Georgescu, B., Kamen, A., Kayvanpour, E., Amr,
A., Sedaghat-Hamedani, F., Haas, J., Katus, H., Meder, B., Hornegger,
J., Comaniciu, D., 2014a. Robust image-based estimation of cardiac tissue
parameters and their uncertainty from noisy data. In: MICCAI. Vol. 8674
of LNCS. Springer, pp. 9–16.
Neumann, D., Mansi, T., Grbic, S., Voigt, I., Georgescu, B., Kayvanpour,
E., Amr, A., Sedaghat-Hamedani, F., Haas, J., Katus, H., Meder, B.,
Hornegger, J., Kamen, A., Comaniciu, D., 2014b. Automatic image-tomodel framework for patient-specific electromechanical modeling of the
heart. In: IEEE Int. Symp. Biomed. Imaging. pp. 935–938.
Neumann, D., Mansi, T., Itu, L., Georgescu, B., Kayvanpour, E., SedaghatHamedani, F., Haas, J., Katus, H., Meder, B., Steidl, S., Hornegger,
J., Comaniciu, D., 2015. Vito – a generic agent for multi-physics model
personalization: Application to heart modeling. In: MICCAI. Vol. 9350
of LNCS. Springer, pp. 442–449.
Nguyen-Tuong, D., Peters, J., 2011. Model learning for robot control: a
survey. Cogn. Process. 12 (4), 319–340.
Noble, D., 2002. Modeling the heart – from genes to cells to the whole organ.
Science 295 (5560), 1678–1682.
Powell, M. J., 2009. The BOBYQA algorithm for bound constrained optimization without derivatives. Cambridge NA Report NA2009/06.
Prakosa, A., Sermesant, M., Allain, P., Villain, N., Rinaldi, C., Rhode,
K., Razavi, R., Delingette, H., Ayache, N., 2013. Cardiac electrophysio30

logical activation pattern estimation from images using a patient-specific
database of synthetic image sequences. IEEE T. Biomed. Eng.
Rosenbrock, H., 1960. An automatic method for finding the greatest or least
value of a function. Comput. J. 3 (3), 175–184.
Schmid, H., Nash, M., Young, A., Hunter, P., 2006. Myocardial material parameter estimation – a comparative study for simple shear.
J. Biomech. Eng. 128 (5), 742–750.
Seegerer, P., Mansi, T., Jolly, M.-P., Neumann, D., Georgescu, B., Kamen,
A., Kayvanpour, E., Amr, A., Sedaghat-Hamedani, F., Haas, J., Katus, H., Meder, B., Comaniciu, D., 2015. Estimation of regional electrical
properties of the heart from 12-lead ECG and images. In: Statistical Atlases and Computational Models of the Heart – Imaging and Modelling
Challenges. Vol. 8896 of LNCS. Springer, pp. 204–212.
Sermesant, M., Billet, F., Chabiniok, R., Mansi, T., Chinchapatnam, P.,
Moireau, P., Peyrat, J.-M., Rhode, K., Ginks, M., Lambiase, P., Arridge,
S., Delingette, H., Sorine, M., Rinaldi, C. A., Chapelle, D., Razavi, R.,
Ayache, N., 2009. Personalised electromechanical model of the heart for
the prediction of the acute effects of cardiac resynchronisation therapy.
In: Functional Imaging and Modeling of the Heart. Vol. 5528 of LNCS.
Springer, pp. 239–248.
Sutton, R. S., Barto, A. G., 1998. Reinforcement learning: An introduction.
Vol. 1. MIT press Cambridge.
Tesauro, G., 1994. Td-gammon, a self-teaching Backgammon program,
achieves master-level play. Neural Comput. 6 (2), 215–219.
Wallman, M., Smith, N. P., Rodriguez, B., 2012. A comparative study of
graph-based, eikonal, and monodomain simulations for the estimation of
cardiac activation times. IEEE T. Biomed. Eng. 59 (6), 1739–1748.
Wallman, M., Smith, N. P., Rodriguez, B., 2014. Computational methods
to reduce uncertainty in the estimation of cardiac conduction properties
from electroanatomical recordings. Med. Image Anal. 18 (1), 228–240.
Wang, V. Y., Lam, H., Ennis, D. B., Cowan, B. R., Young, A. A., Nash,
M. P., 2009. Modelling passive diastolic mechanics with quantitative mri
of cardiac structure and function. Med. Image. Anal. 13 (5), 773–784.
Wang, Y., Georgescu, B., Chen, T., Wu, W., Wang, P., Lu, X., Ionasec, R.,
Zheng, Y., Comaniciu, D., 2013. Learning-based detection and tracking
in medical imaging: a probabilistic approach. In: Deformation Models.
Vol. 7 of LNCVB. Springer, pp. 209–235.
31

Westerhof, N., Elzinga, G., Sipkema, P., 1971. An artificial arterial system
for pumping hearts. J. Appl. Physiol. 31 (5), 776–781.
Wong, K. C., Sermesant, M., Rhode, K., Ginks, M., Rinaldi, C. A.,
Razavi, R., Delingette, H., Ayache, N., 2015. Velocity-based cardiac contractility personalization from images using derivative-free optimization.
J. Mech. Behav. Biomed. 43, 35–52.
Xi, J., Lamata, P., Niederer, S., Land, S., Shi, W., Zhuang, X., Ourselin,
S., Duckett, S. G., Shetty, A. K., Rinaldi, C. A., Rueckert, D., Razavi,
R., Smith, N. P., 2013. The estimation of patient-specific cardiac diastolic
functions from clinical measurements. Med. Image. Anal. 17 (2), 133–146.
Zettinig, O., Mansi, T., Neumann, D., Georgescu, B., Rapaka, S., Seegerer,
P., Kayvanpour, E., Sedaghat-Hamedani, F., Amr, A., Haas, J., Steen,
H., Katus, H., Meder, B., Navab, N., Kamen, A., Comaniciu, D., 2014.
Data-driven estimation of cardiac electrical diffusivity from 12-lead ECG
signals. Med. Image Anal., 1361–1376.
Zheng, Y., Barbu, A., Georgescu, B., Scheuering, M., Comaniciu, D., 2008.
Four-chamber heart modeling and automatic segmentation for 3-D cardiac
CT volumes using marginal space learning and steerable features. IEEE
T. Med. Imaging 27 (11), 1668–1681.

32

