Occupation measure methods for modelling
and analysis of biological hybrid systems
Alexandre Rocca ∗,∗∗ Marcelo Forets ∗ Victor Magron ∗
Eric Fanchon ∗∗ Thao Dang ∗
∗

arXiv:1710.03158v2 [] 15 Mar 2018

Univ. Grenoble Alpes, CNRS, Grenoble INP, VERIMAG, 700,
avenue centrale,38000, Grenoble, France
∗∗
Univ. Grenoble Alpes, Grenoble 1/CNRS, TIMC-IMAG, UMR 5525,
38041, Grenoble,France 1
Abstract: Mechanistic models in biology often involve numerous parameters about which we
do not have direct experimental information. The traditional approach is to fit these parameters
using extensive numerical simulations (e.g. by the Monte-Carlo method), and eventually revising
the model if the predictions do not correspond to the actual measurements. In this work we
propose a methodology for hybrid system model revision, when new types of functions are needed
to capture time varying parameters. To this end, we formulate a hybrid optimal control problem
with intermediate points as successive infinite-dimensional linear programs (LP) on occupation
measures. Then, these infinite-dimensional LPs are solved using a hierarchy of semidefinite
relaxations. The whole procedure is exposed on a recent model for haemoglobin production in
erythrocytes.
Keywords: biological modelling, hybrid dynamical system, optimal control problem,
semidefinite optimization, occupation measures.
1. INTRODUCTION AND CONTEXT
Mechanistic models in biology generally involve many parameters. The value of a given parameter can be either
measured directly in a dedicated experiment (e.g. measurement of a kinetic parameter of a biochemical reaction
in enzymology), or inferred from data which provide relationships between parameters and other biological entities.
In this paper, we work with biological mechanisms modelled by ordinary differential equations (ODEs) and hybrid
dynamical systems. The motivation for using such models
is to give quantitative predictions, when sufficient data is
available to validate the model. And whenever new data
and knowledge of various types become available, they can
be incorporated within a formal framework.
A basic issue in biological systems modelling is the determination of numerical values for the parameters, or more
generally a subset of the parameter space, under which the
model agrees to some extent with the available data. We
focus on multiple-step experiments, in which a biological
system is perturbed or measured during its evolution.
In the biological systems modelling literature, it is common to synthesise parameters using a Monte-Carlo sampling of the parameter space, which is validated then by
numerous simulations. An important effort to formalize
and validate the parameter synthesis of biological models
has been made in works such as Donzé (2010); Mobilia
(2015); Dreossi (2016); Beneš et al. (2016), or Rumschinski
et al. (2010). Other such as Cardelli et al. (2017) or
1 For all authors emails follow forename.surname@univ-grenoblealpes.fr

Bartocci et al. (2013) design ODE models satisfying sets
of temporal constraints. When model simulation does not
reproduce satisfactorily available experimental data, to a
degree which depends on data quality, for any admissible
parameter value, the model has to be revised. One way
of revising the model is to represent parameters using
different types of functions of time, reflecting underlying
biological mechanisms. We introduce a systematic way,
based on formal methods, to study mechanistic biological models in their experimental context and revise parameters to produce conservative results with respect to
experimental data. In this work, we consider a problem
of model revision, defined as finding time varying laws of
parameter evolution that minimizes the error in matching
experimental measurements. Concretely, it is the following
optimization problem:
nexp

inf
(x, u)

X

dist(m(x(Tj )), zj )

(1)

j=1

where x is a vector of biological variables, such as concentrations, whose evolution is modelled by trajectories of a
biological dynamical system. Time varying parameters are
represented by the input variables u (modelling biological
parameters) such that ∀t ∈ [0, T ] , u(t) ∈ U. X0 is the set
of initial conditions and the set of pairs {(Tj , zj )}j is the
set of data points, for 1 ≤ j ≤ nexp , in the time frame
[0, T ]. An experimental measurement is a function of the
variables x and is modelled via the function m(x).
The framework of our approach is a mathematical formalization of experimental protocols as hybrid dynamical systems, describing biological systems of interest and
experiments which are performed on them. However, the

algorithm we provide can be applied to any biological
hybrid system with similar model revision problems.
In this paper we address the parameter synthesis and
model revision problems (1) by formulating a particular
instance of the optimal control problem with intermediate
costs, when the objective function depends on the systems
trajectory and control inputs at a given set of time
points. This problem is then approximated by multiple
hybrid optimal control problems (HOCP) with one final
cost. To this end, we apply a recently developed method
of Zhao et al. (2017) from the field of certified convex
optimization to globally solve these HOCP. However, The
method described in Zhao et al. (2017) produces piecewise
optimal control functions which either may not correspond
to biological knowledge of parameter variations or may
be difficult to yield coherent and meaningful biological
interpretations. Consequently, in order to respect realistic
constraints on parameters, we use smooth approximations
of the generated control input, in order to revise the given
model while maintaining good data fitting accuracy.
The method is demonstrated on a hybrid system modelling
haemoglobin production presented in Section 4.1.
1.1 Related work
The hybrid formalism has previously been used as an abstraction method to simplify complex mechanisms which
are hard to analyse as seen in Noel et al. (2011); Rocca
et al. (2016), or to represent “jump” evolution such as activation processes in genes regulatory networks for example
using the stochastic formalism as in Li et al. (2017).
Optimal control theory and variation theory have been
applied to biological systems in several works. Most of
them address the classical problem of finding a correct
input such that the system reaches a desired state. For
example, one can control drug input such that a patient
attains a healthy state, see Ledzewicz and Schättler (2007)
or Caraguel et al. (2016). Another example is the control
of some input in population studies as detailed in Bodine
et al. (2008). A detailed review on the use of optimal
control in systems biology can be found in Lenhart and
Workman (2007). The problem of parameter estimation
in presence of multiple data, also called data assimilation,
is stated in (Lenhart and Workman, 2007, Chapter 26).
The optimal control problem for specific classes of hybrid systems has been investigated in several domains,
such as mechanical systems in Pace and Burden (2017)
and switched-mode systems in Wardi et al. (2015); Xu
and Antsaklis (2004); Bengea and DeCarlo (2005). More
generally, Pakniyat and Caines (2014) relies on Dynamic
Programming and an extension of Pontryagin’s Maximum
Principle. However, these approaches need a-priori knowledge either on the sequence of discrete transitions, or on
the number of visited subsystems. To perform optimal
control on hybrid systems, we build our work on the techniques from Zhao et al. (2017), which proposes a method
to obtain a global solution for hybrid systems with statedependent transitions, without any a-priori knowledge on
the execution and the sequence of transitions. We refer to
(Zhao et al., 2017, Section 1.1) and references therein for
more details on optimal control of hybrid systems.

Semidefinite programming (SDP) eases the resolution of
hard optimization problems and yields conservative results ensured by positivity certificates. In Lasserre (2001),
hierarchies of semidefinite relaxations were introduced
for static polynomial optimization. The definition of an
infinite-dimensional linear program (LP) over occupation
measures, for optimal control problems, was first introduced in Vinter (1993). From this infinite-dimensional LP,
Lasserre et al. (2008) defines hierarchies of Linear Matrix
Inequalities (LMI) relaxations, to synthesise a sequence
of polynomial controls converging to the solutions of the
optimal control problem. In Abdalmoaty et al. (2013) the
authors propose an extension to piecewise affine systems.
Our underlying idea of constructing a suboptimal control
with an iterative algorithm is similar to (Abdalmoaty
et al., 2013, Section 4). However, we use this scheme to
find input functions allowing to reproduce data not only
at a final time point but also at intermediate time points.
We make use of the recent method proposed in Zhao
et al. (2017), which relies on occupation measures and a
sequence semidefinite relaxations to produce a sequence of
polynomial controls converging to the optimal solution of
a HOCP. There exist other methods which use occupation
measures and LMI relaxations to produce both admissible
controls and converging outer-approximations of the backward reachable set (BRS) Majumdar et al. (2014); Shia
et al. (2014), or the region of attraction (ROA) Korda
et al. (2014).Finally, we note that finding a sequence of
converging outer-approximations for all valid parameters
sets, such as in Streif et al. (2013); Dreossi (2016), is
another crucial issue in the context of systems biology.
When dealing with hybrid systems, Shia et al. (2014) can
be applied to solve this problem.
1.2 Main contributions
Given a hybrid system modelling a multiple-step biological protocol and a set of experimental measurements,
we propose a numerical algorithm to solve the model
revision problem. This algorithm generates an admissible
input function such that the revised model reproduces the
experimental measurements, and can be interpreted. We
release a software package 2 , in MATLAB, implementing
this algorithm. We evaluate this algorithm and its implementation, on the haemoglobin production model taken
from Bouchnita et al. (2016).
The rest of the paper is organized as follows: in Section 2
we give the necessary notations on dynamical hybrid systems and the optimal control problem. Section 3 presents
our main contribution to the resolution of the hybrid
optimal control problem with intermediate points. Finally,
in Section 4 we study and discuss the model revision of
the haemoglobin production model taken from Bouchnita
et al. (2016).
2. PRELIMINARIES
We first give the notations and basic notions on controlled
hybrid systems, as well as the definition of the HOCP.
Given x ∈ Rn , let xi denote its i-th component. In general,
letters in bold font denote multidimensional elements, and
2

https://github.com/biosdp/biosdp

normal font unidimensional ones. Let B := {true, false}
be the set of Booleans. Let R[x] denote the ring of real
polynomials in x ∈ Rn , and let Rd [x] be the subspace of
polynomials whose degree is at most d. Let T be the time
interval [0, T ], where T is the final time (possibly ∞).

Given a CHS H, the hybrid optimal control problem is
defined as follows. Let X0 , and XT , be the initial set, and
target set defined by:
a
a
X0 :=
X0,i , XT :=
XT,i ,
(4)

Consider the n-dimensional ODE with inputs:
ẋ(t) = f (t, x(t), u(t)),
(2)
n
m
n
with f : T ×R ×R → R a vector field which is Lipschitz
continuous in x and piecewise continuous in u. Let X and
U be compact subsets of Rn and Rm respectively. Here,
u : T → U is a feasible input function which represents
time varying parameters, or external commands. The tuple
F := (T , X, U, f )
(3)
defines a continuous dynamical system.

where X0i and XT i are a compact subsets of Xi for
each mode i ∈ I. Let i0 and iT be the initial mode
and the final mode at time T , respectively. Then, given
(i0 , x(0)) = (i0 , x0 ) ∈ X0 and u : T → U an input
function, we say that for T > 0, (x(t), u(t)) ∈ P is an
admissible pair on T and P is the set of admissible pairs,
if (i, x(t)) ∈ X is a trajectory of CHS as defined by (Zhao
et al., 2017, Algorithm 1) and (iT , x(T )) ∈ XT . Finally,
the hybrid optimal control problem for a CHS H, (HOCP),
is defined by:
Z T

∗
Jhocp
hλ(x(t)) t, xλ(x(t)) (t), u(t) dt
:=
inf
(5)
(x,u)∈P 0

+ Hλ(x(T )) xλ(x(T )) (T ) .

We recall the definition of controlled hybrid systems a
dynamic hybrid systems formalism defined in Zhao et al.
(2017).
Definition 1. (Controlled hybrid system). A controlled hybrid system (CHS) is defined by the tuple: H =
(I, E, D, U, F, S, R) where:
• I ⊂ N is the set of mode indices, and nmodes the
number of modes.
• E ⊆ I × I is the set of transitions e = (i, j) between
two modes: i is the source mode, and j the destination
mode.`
• X := i∈I Xi is the disjoint union of domains of H
and Xi the domain of the mode i. We note that Xi
the dimension
is a compact subset of Rni with ni `
of the mode i. The disjoint union
can simply
be considered as a labelling operation on the set of
domains by I, that is the set of mode indices.
• U is the set of input values of H.
• F := {Fi }i∈I is the set of continuous dynamical
sub-systems associated to each mode. The dynamical
system associated to mode i is:
Fi := (T , Xi , U, fi ) ,
with fi : T × Xi × U → Rni a vector field polynomial
in x and
` affine in u.
• S := e∈E Se is the disjoint union of guards Se ⊆ Xi
associated to each transition e = (i, j) ∈ E. The guard
S(i,j) defines the switch condition from i to j: for
x ∈ Xi , if x ∈ S(i,j) then the system at x can make
the transition from mode i to mode j.
• R := {Re }e∈E is the set of reset maps, each reset
map Re : Se → Xj being associated to a transition
e := (i, j) ∈ E and it defines how the continuous
variables may change after the discrete transition
from mode i to mode j.
Additionally, the CHS defined in Zhao et al. (2017) must
respect a few assumptions. All the guards S(i,·) are disjoint,
and S(i,j) ⊆ ∂Xi , with ∂Xi designating the border of Xi ,
for each pair of modes i and j. The initial set is restricted
to a single point x0 , with an associated mode i0 . The
vector fields fi are affine in u and have a nonzero normal
component on the boundary of Xi . These assumptions
ensure that any CHS is deterministic. Noting λ(x(t)) the
function which associates to an instantaneous state x(t)
its corresponding mode, these assumptions ensure that the
mode corresponding to x(t) is unique.

i∈I

i∈I

where {hi : [0, T ] × Rni × Rm → R}i∈I and {Hi : Rni →
R}i∈I are set of measurable functions.
3. OPTIMAL CONTROL FOR MODEL REVISION
In our biological context hybrid dynamical systems can
model multiple steps experiments. Indeed, experimental
protocols associated to these experiment can be considered
as a set of concurrent processes where each process is
modelled as a hybrid system. Thus, the hybrid system of
the protocol can be described as the parallel composition
of the hybrid system describing each process. It is then
crucial to proceed to model revision while taking into
account the biological system in the evolving environment
of the complete protocol. For this purpose, we propose in
this section a method to produce time varying parameters
reproducing multiple data in the context of a multiplephase protocol modelled by a hybrid system. However, we
also argue that this method can be use for more general
biological systems modelled as CHS. We solve the model
revision problem of a dynamical hybrid system modelling
a biological system together with a set of experiments.
Therefore, we search for parameters as time varying functions fitting a set of data points, defined in the introduction
as in (1). In this aim, Section 3.1 first formulates (1) as
a particular instance of the optimal control problem on
hybrid systems with intermediate costs. Then, we propose
a first approximation as a set of instances of the HOCP
(5) defined previously. The solution of each sub problem is
obtained using the previous results from (Zhao et al., 2017,
Section 4). Finally, in Section 3.2 we explain the complete
algorithm addressing our initial problem.
3.1 Problem statement
We provide a method to find time varying parameters
of biological system, modelled as input functions u(t), in
order to fit the hybrid systems model to a set of experimental data. Thus, we write our problem as an optimal
control problem where desired input functions are the
optimal controls which minimize the distance of the results
produced by the model and these experimental data. Experimental measurements, represented by a function m(x),

are performed at given specific times Tj , 1 ≤ j ≤ nexp .
Let zj be the observed value of the experimental measurement at time Tj , then nexp is the number of experimental
data points.
` Let XTj,i be compact subsets of Xi , and
XTj := i∈I XTj ,i . As in (5), let (i0 , x(0)) ∈ X0 , and
suppose that we are given a set of time values {Tj }, with
1 ≤ j ≤ nexp , and Tnexp = T . Given an input function
u : T → U, we say that (x, u) ∈ Pint is an admissible
pair for a problem with intermediate points and Pint the
associated set of admissible pairs, if (i(t), x(t)) ∈ X is a
CHS trajectory, and (iTj , x(Tj )) ∈ XTj for all j.
Let H(x(Tj )) be a cost at time Tj , and h(t, x(t), u(t))
a running cost for the whole [0, T ] interval. The optimal
control problem with intermediate points for the CHS H
is then:
Z T
∗
h(t, x(t), u(t))dt
Jhocp
:=
inf
(x,u)∈Pint 0
(6)
X
H (x(Tj ))
+
0≤j≤nexp

In our biological context, H (x(Tj )) = ||m(x(Tj )) − zj ||22 .
To our best knowledge there is no method to efficiently address directly problem (6), and obtain converging sequence
of solution in a similar manner to the simpler problem
(5) address in Zhao et al. (2017). Thus, we search for
an admissible solution using a greedy algorithm. Consequently, this solution is a good trade-off between computation cost, optimality and flexibility as we can handle
a large class of models. Moreover, this method does not
constraint the form of the sought times parameters: this
reduces the assumptions during the modelling and eases
the final biological interpretations, unlike methods based
on simulations. Given 1 ≤ j ≤ nexp , let
Z Tj
Jj (t, x(t), u(t)) :=
h(t, x(t), u(t)) + H (x(Tj )) ,

mode i for the next optimization on the interval [Tj , Tj+1 ].
Let ũ(t) and (i(t), x(t)) be respectively the control and
the trajectory, for t ∈ [0, T ]. They are respectively defined
by the concatenation of all the controls ũ(j) (t) and the
trajectories (i(j) (t), x(j) (t)) on the sub-intervals [Tj−1 , Tj ].
By construction, (x(t), ũ(t)) is an admissible pair for (6),
(j)
as (iTj , x(Tj )) = (iTj , x(j) (Tj )) ∈ XTj .
Remark 1. We emphasize that (x(t), ũ(t)) is not necessary
an optimal solution for (6). Moreover, as the optimization
problem (7) is obtained through a greedy scheme, we have
no guarantee that its optimal cost Jj∗ is inferior to a given
ε. However, as we search to equally fit all the data points
searching iteratively for the control is satisfiable solution.
3.2 Implementation
Let (Tj , zj ), 0 ≤ j ≤ nexp be pairs of experimental data
points and their measurement time, and we also note
i0 , and x0 the initial mode and initial conditions of the
studied CHS H respectively. Let r be a given starting
relaxation degree. Algorithm 1 finds an admissible solution
to (6), by solving a sequence of the HOCP (7) for each
experimental data point (Tj , zj ). For each j, the degree of
the polynomial control ũ(j) (x(t), t) is determined as the
smallest degree such that ||m(x(Tj )) − zj ||22 ≤ ε. Indeed,
in the context of biological system modelling we desire to
obtain a control of degree as small as possible to avoid
overfitting. We note that as explained in Remark 1, we
cannot ensure the converge to a solution of accuracy ε.
Thus, ε is only a stopping criteria. Then, for each iteration
over j, Algorithm 1 is decomposed in three steps.

The first step is the procedure HOCP, associated to an
instance of the HOCP (5) for j-th pairs (Tj , zj ). Given
a relaxation order dr ≥ r, we solve the relaxed primal
defined in (Zhao et al., 2017, Section 5.1). From Zhao et al.
Tj−1
(2017) method we first obtain Mdr (yµi ), the sequence mowith T0 = 0, and Tnexp = T , such that
ment matrices of degree dr , associated to the occupation
X
(d )
J(t, x(t), u(t)) =
Jj (t, x(t), u(t)).
measure µi of each mode i ∈ I. We also obtain J j r an
1≤j≤nexp
under approximation of the optimum of (7). The second
(j)
(j)
Noting (i (t), x (t)) a trajectory of a CHS H on the step is the(j)procedure Synth, which returns the admissible
interval Tj := [Tj−1 , Tj ], and similarly ũ(j) (t) the control control ũ (t, x) of degree du ≤ dr using a truncated
on Tj , we consider the following problem as particular moment matrix Mdu (yµi ) of Mdr (yµi ) at the reduced
degree du . The third and last step is the procedure Simu.
instance of (5):
It performs the validation that the synthesised control ũ(j)
Jj∗ :=
inf Jj (t, x(j) (t), ũ(j) (t))
yields ||m(x(Tj ))−zj ||22 ≤ ε. This step is done by approxi(x(j) ,ũ(j) )
mating the trajectory of the controlled hybrid system using
s.t.
a solver of ODE with discrete events to produce numerical
simulations. If in iteration j, ||m(x(Tj )) − zj ||22 ≤ ε, then
x(j) cont. part of a trajectory (i(j) , x(j) ) on Tj ,
x(j) (Tj ) and the corresponding mode if reached at t = Tj
ũ(j) (t) ∈ U , ∀t ∈ Tj ,
by the numerical simulations are the initial conditions for
(i(j) (t), x(j) (t)) ∈ X , ∀t ∈ Tj ,
the next iteration j + 1. Otherwise, the Ctrl Synth and
if j = 1 ,
Simulate procedures are repeated while increasing the
degree of the synthesised polynomial control until du = dr .
(i(1) (0), x(1) (0)) ∈ X0 ,
In case the condition ||m(x(Tj )) − zj ||22 ≤ ε is still not
if j ≥ 2.
satisfied, the relaxation order dr is increased, and the three
(i(j) (Tj−1 ), x(j) (Tj−1 )) = (i(j−1) (Tj−1 ), x(j−1) (Tj−1steps
)) , are repeated.
(i(j) (Tj ), x(j) (Tj )) ∈ XTj .

(d )

(7)
We note that if a transition i → i0 occurs at the time Tj of
the interval [Tj−1 , Tj ], we retain only the left part in the

If ε ≤ J j r then we are sure that we the given initial condition at step j, there is no control such that ||m(x(Tj )) −
zj ||22 ≤ ε. Consequently, we keep our previous result ũ(j)
and the corresponding mode if reached at t = Tj by the

numerical simulations are the initial conditions for the
next iteration j + 1.
Algorithm 1.
1: procedure Algorithm 1(H, {(Tj , zj )}j , i0 , x0 , ε, r)
2:
Tinit = 0
3:
for all experimental data (Tj , zj ) do
4:
du = 0, dr = r, err = +∞
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

(dr )

while err ≥ ε ∧ J j

≤ ε do

(dr )

Jj

, Mdr (yµ ) = HOCP(H, i0 ,...
...x0 , Tinit , Tj , zj , dr )
while err ≥ ε and du ≤ dr do
ũ(j) (x(t), t) = Synth(Mdr (yµ ), du )
(if , x(j) (t)) = Simu(H, ũ(j) (x(t), t),...
... i0 , x0 , Tinit , Tj )
err = H(x(j) (Tj ), zj )
increase du
end while
increase dr
end while
i0 = if
x0 = x(j) (Tj )
Tinit = Tj
end for
end procedure

4. RESULTS AND DISCUSSION
In this section, using the method developed in Section 3,
we revise the model of haemoglobin production by finding
a better fit for the time varying parameter noted k3 in
Bouchnita et al. (2016), with respect to the same error
function. To be consistent with our notations from Sections
2 and 3, we will note u(t) the input modelling the time
varying parameter k3 .
4.1 Application to haemoglobin production and results
The paper Bouchnita et al. (2016) addressed, on an experimental protocol model, the problem of refining the parameters space, and to fit multiple data sets of experimental results. In Bouchnita et al. (2016), the experimental
protocol was not explicitly formalized as a hybrid system.
However, we use the parameters value from this previous
work as starting instantiation to provide less restrictive
constraints on the time varying parameter k3 .
The ODEs (fctrl ) model the evolution of the haemoglobin
production in the differentiating erythrocyte cells situated
in the bone marrow. The variables x1 to x4 represent
respectively the internal iron in the cell Fe, the heme
H , the globin G, and the haemoglobin Hb. The hybrid
system H models an experimental protocol designed to
measure the integration of iron inside heme (H) at multiple
milestones of the cell differentiation. For example the data
point at time t = 7 hours in Table 3, is obtained through
the following procedure: we first start with a control batch
of cells, then at time t = 4 hours after the start of
the experiment, the culture medium is perturbed with an
injection of measurable radioactive iron 59Fe for a subset
of the cells. This perturbation implies new ODEs (frad )

Table 1. Dimensions (with xc ), vector fields,
domains, and input sets for the controlled hybrid system H of the haemoglobin production
model.
Mode

ni

i=1
i=2

fi (t, x, u)

5

Xi

fctrl (t, x, u)

9

frad (t, x, u)

Ui

[0, 4] ×

[0, 1]4

[0, 1]

[4, 7] ×

[0, 1]8

[0, 1]

[7, 8] ×

[0, 1]4

[0, 1]

i=3

5

fctrl (t, x, u)

.
i = ..

..
.

..
.

i = 13

5

fctrl (t, x, u)

[45, 52] × [0, 1]4

[0, 1]

frad (t, x, u)

[0, 1]8

[0, 1]

i = 14

9

..
.

..
.

[52, 55] ×

Mode

e = (i, j)

Se

i=1

(1, 2)

t == 4

i=2

(2, 3)

t == 7

i=3

(3, 4)

t == 8

.
i = ..

..
.

..
.

i = 13

(13, 14)

t == 52

Re

h

I5,5
O4,4

i



I5,5 , O4,4

h

I5,5
O4,4

i

I5,5
O4,4

i



..
.

h

Table 2. Transitions, guards, and reset maps of
the controlled hybrid system H.
modelling the evolution of two inter-dependant models: the
model of non-radioactive haemoglobin production and the
model of haemoglobin production with radioactive species.
Three hours after the perturbation with radioactive iron,
the total radioactive heme is measured, meaning the heme
free in the cell and the one in the radioactive haemoglobin.
This measurement is given by the formulas 59H + 4 59Hb.
The data in Table 3 is then the observed radioactivity divided by three hours. Finally, these measurements provide
results on the variation during the cell differentiation of
the integration of iron in heme, which is associated to the
parameter k3 .
The controlled hybrid system H associated to the experiment studied in Bouchnita et al. (2016) and the
haemoglobin production model are given, in a shortened
version, in Tables 1 and 2. The ODEs (fctrl ) and (frad ), as
well as the numerical values of the parameters, are given
in Appendix. In the implementation, we also introduce a
variable xc modelling time, whose derivative is equal to 1.
For numerical reasons, it is necessary to scale the parameters and state variables, making it easier for the solver
to succeed in solving the relaxed problem. Similarly, to
facilitate the numerical optimization we rewrite the control
variable u(t) ∈ U = [0, 1] as u(t) = ζ û(t), with ζ  1 and
û(t) ∈ [0, 1/ζ]. While the scale factor ζ may take different
values depending on the numerical optimization details,
the objective control u(t) always evolves in [0, 1].
Now that we have a valid controlled hybrid system, H, we
solve the optimal control problem with intermediate time
points defined in (6), using the method from Section 3.1

2.5

×10-4

vein, we add another penalization cost c2i (t) = (u(Tj ) −
u(t))2 to avoid when the control strongly varies between
two iterations j on the interval [Tj−1 , Tj ] and j + 1 on
[Tj , Tj+1 ] (with the exception of the first iteration). This
leads to hi (t, x(t), u(t)) = c1i (t) + c2i (t). Let us note that,
even if these additional costs can eventually degrade the
accuracy of the data fitting, we gain in terms of biological
interpretation of the resulting traces.

Generated optimal control and associated approximating functions

2.25

2

1.75

k 3(t) [hours-1 ]

1.5

1.25

1

0.75

0.5
Generated Optimal control u gen (t)
Piecewise polynomial u poly (t)

0.25

Hill function u hill(t)
Step function u step (t)
0
0

5

10

15

20

25

30

35

40

45

50

55

60

Time [hours]

Fig. 1. Synthesised optimal control and various approximations that yield a realistic interpretation.
Fe59(t)

0.08

Gtot(t)

0.12

0.06

Gtot(t) [mol/m 3 ]

Fe59(t) [mol/m3 ]

0.1

0.04

0.02

0.08
0.06
0.04
0.02

0

0
0

10

20

30

40

50

60

0

10

20

H59(t)+4*Hb59(t) [cpm/(1e-7 L)h-1 ]

Time [hours]

30

40

50

60

Time [hours]

Observed variable: H59(t)+4*Hb59(t)

500
400
300
200
100
0
0

10

20

30

40

50

60

Time [hours]

Fig. 2. Radio-active variables 59Fe, Gtot corresponding to
x5 and x7 in (frad ), as well as, the comparison of the
measurement function results to the data.
and its implementation in Section 3.2. The experimental
measurement is modelled by the function m(x) = x6 +4x8 .
Thus, we set H(x(Tj )) := (x6 (Tj ) + 4x8 (Tj ) − zj )2 , as we
search to minimize the total residual error term:
p
X
H(x(Tj ))
P
.
(8)
εtotal =
1≤j≤nexp zj
1≤j≤nexp

The original experimental data points (Tj , zj ) are given in
Table 3. Here, the input control u(t) models some hidden
Time (h)
cpm
Measure ( 1e−7
)
L·h−1

7
16

11
85

19
348

27
391

35
399

45
481

55
395

Table 3. Experimental data points (Tj , zj ) used
as references.

Finally, by partitioning the computation in the time domain, we can greatly reduce the computational cost at
each iteration. More technically, since the transitions of
the hybrid system H are fully determined by the time
t, we can pre-compute the function λ : R+ → I, which
associates a mode λ(t) to each time instant t. Thus, each
iteration j of Algorithm 1 can be restrained to the subhybrid system Hj of H, constituted by the modes visited
in the interval [Tj−1 , Tj ].
For numerical implementation, the problem on measures
is formulated in SPOTLESS 3 , and then we extract the
primal solution provided by a primal-dual SDP solver. To
do so, we use the implementation from Zhao et al. (2017)
to generate the dual problem. As an SDP solver we used
MOSEK Andersen and Andersen (2000) v.7.1. These tools
are used in MATLAB v.9.0 (R2016a). Performance results
are obtained with an Intel Core i7-5600U CPU (2.60 GHz)
with 16Gb of RAM running on Debian 8. We only solve
the problem for a relaxation order r = 4, as any higher
order would be too memory expensive. In this particular
example as constant input leaded to satisfying results, we
did not impose any constraint  in the algorithm. Using
this configuration, the total time taken by Algorithm 1
is 2107s, with 1700s spent in the HOCP procedure, and
390s in the Synth procedure. On Figure 1, the control
generated by Algorithm 1 is shown in blue. This control
is piecewise, and clearly divided in two phases: before and
after t equals 11 hours. However, the control synthesised
is still difficult to interpret as a biological phenomenon.
Consequently, we propose three additional fits of this
control to ease interpretation by using functions closer to
biological knowledge. In Table 4 one can find the total
error associated to all the possible controls, as well as the
previous result of Bouchnita et al. (2016). In Figure 2,
we show a graphical representation of how closely each
function can control the model to reach the desired data
points.

mechanism which evolves with the differentiation of the
cells. It should be the same function of time for both the
control and the radioactive cells batch. However, as the
control generated by Algorithm 1 is piecewise for each
mode, and the fact that our data are on the radioactive
species only, the solution of the optimization problem with
only a final cost H(x(Ti )) is not balanced, having a much
stronger control in the modes where the radioactive species
are evolving.
A workaround for the balancing problem is the following.
We add a small penalization cost c1i (t) = (0.01u(t))2 to
equilibrate the control when i corresponds to a mode
with radioactive species, otherwise c1i (t) = 0. In a similar

Control Type

εtotal

Original Control Bouchnita et al. (2016)
Algo 1 generated
Step function fit
Piecewise Polynomial fit
Hill function fit

0.23
0.096
0.12
0.13
0.075

Table 4. Total error εtotal associated to each
proposed input.

3

The
SPOTLESS
implementation
https://github.com/spot-toolbox/spotless

was

taken

from

4.2 Discussion
In a simulation-based approach, we have to propose for
the desired time varying parameter, a template function
to fit the data, e.g. a polynomial of given degree. If we
want to fit a polynomial of higher degree, the simulations
have to be run again multiple times. On the contrary, the
proposed approach returns a control signal, and since the
fit to data points is performed a posteriori, there is no
additional computation cost in refining the model.
From the form of the experimental data points, an usual
hypothesis is that u(t) should be similar to a jump function, with a low value for the two first points, and a
higher one for the following ones. However, even with
such information a good fit is not easily achieved with
simulations.
The control generated with Algorithm 1 returns the expected “jump” behaviour for u(t), and even with a low
relaxation degree, the total residual error for the generated
control is 9.59% which is much lower than the 22.8% from
Bouchnita et al. (2016).
We first fit a step function to the generated control, with
a change at t = 11. The associated error of 12.24% is
still lower than Bouchnita et al. (2016), yet being higher
than the generated control mainly due to the second-tolast point.
The second fit is a piecewise polynomial function in two
pieces. The first piece, for t ∈ [0, 11], is a polynomial of
degree 2 while the degree of the second, for t ∈ [11, 55], is
4. This proposed input control allows to reproduce more
accurately, than the step function, the third data point.
However, its accuracy is worse on the first and two last
points. The total error associated to this control is 13%,
being overall the worst of the proposed fits.
Lastly, we fit a Hill function, a function used to model
the kinetics of a class of biochemical reactions and which
is a very common way to represent biological activation
processes. The associated total error is 7.5%, which is
the lowest, taking advantages from both the step function
and the piecewise polynomial function. In this case, the
inaccuracy also mainly comes from the second-to-last
point, which is quite separated from trend of the other
experimental points, and may be due to some experimental
problems (no standard deviation results were available).
Without taking this point into consideration for the error
computation the error falls to 3% for the Hill function fit.
On this particular example, this method provided a way to
generate a control satisfying intermediate points without
any a priori on a particular form, avoiding the need for
extensive numerical simulations. The generated control is
accurate, and computed in a reasonable time (∼35min),
even for a large hybrid system of 14 modes with at
most 9 continuous variables. Using some fitting functions
afterwards, it is even possible to refine the results and
obtain a more interpretable function for the desired time
varying parameters.
As, in this model, the sequence of transitions is known in
advance, the use of Zhao et al. (2017) to solve (5) at each
iteration of Algorithm 1 is arguable, as other methods can
handle this problem. If needed Algorithm 1 can easily be

adapted with another method to solve the optimal control
problem on hybrid system, using for example simulations.
However, Algorithm 1 in its current form does not require
any knowledge on the sequence of transition and can be
applied to a larger set of biological models.
5. CONCLUSION
In this work, we address an important problem arising in
biological modelling: model revision. We propose a method
for revising an experiment modelled by a hybrid system,
given a set of experimental data points. The method scales
even on large hybrid systems such as the haemoglobin production model, while providing an accurate result, and a
meaningful interpretation, as an activation process, for the
mechanism underlying the revised parameter. The CHS
formalism is motivated by the development of an automatic, and formal modelling of multiple-step experimental
protocols, and to develop new methods for their analysis.
Such formal representations had already been used as
alternative, non-ambiguous languages, in contrast with the
natural language, for the description of experiments Soldatova et al. (2008). However, those works do not consider an
underlying mechanistic model in the form of ODEs. In a
future work, we plan to provide an automatic and rigorous
way to obtain a hybrid system modelling an experimental
protocol with their underlying mechanisms. Moreover, we
plane to investigate, using the CHS formalism, two other
relevant problems in biological systems modelling: finding valid subsets of the parameters space fitting multiple
data points (as an extension Shia et al. (2014)), and the
validation of biological experiments. The focus of this
work has been on ODEs, but multicellular systems and
transport processes are described by partial differential
equation (PDE) models. The extension of semidefinite
programming techniques to PDEs Mevissen et al. (2011)
and their application to biological models would require
further theoretical and numerical developments.
REFERENCES
Abdalmoaty, M.R., Henrion, D., and Rodrigues, L. (2013).
Measures and lmis for optimal control of piecewise-affine
systems. In ECC.
Andersen, E.D. and Andersen, K.D. (2000). The Mosek
Interior Point Optimizer for Linear Programming: An
Implementation of the Homogeneous Algorithm.
Bartocci, E., Bortolussi, L., and Nenzi, L. (2013). A
temporal logic approach to modular design of synthetic
biological circuits. In CMSB.
Beneš, N., Brim, L., Demko, M., Pastva, S., and Šafránek,
D. (2016). Parallel smt-based parameter synthesis with
application to piecewise multi-affine systems. In ATVA.
Bengea, S.C. and DeCarlo, R.A. (2005). Optimal control
of switching systems. automatica.
Bodine, E.N., Gross, L.J., and Lenhart, S. (2008). Optimal
control applied to a model for species augmentation.
Mathematical biosciences and engineering.
Bouchnita, A., Rocca, A., Fanchon, E., Koury, M., Moulis,
J., and Volpert, V. (2016). Multi-scale modelling of
erythropoiesis and hemoglobin production. JIOPM.
Caraguel, F., Lesart, A.C., Estève, F., van der Sanden,
B., and Stéphanou, A. (2016). Towards the design of
a patient-specific virtual tumour. Computational and
mathematical methods in medicine.

Cardelli, L., Češka, M., Fränzle, M., Kwiatkowska, M.,
Laurenti, L., Paoletti, N., and Whitby, M. (2017).
Syntax-guided optimal synthesis for chemical reaction
networks. In CAV.
Donzé, A. (2010). Breach, a toolbox for verification and
parameter synthesis of hybrid systems. In CAV.
Dreossi, T. (2016). Sapo: Reachability computation and
parameter synthesis of polynomial dynamical systems.
arXiv.
Korda, M., Henrion, D., and Jones, C.N. (2014). Controller
design and region of attraction estimation for nonlinear
dynamical systems. IFAC.
Lasserre, J.B. (2001). Global optimization with polynomials and the problem of moments. SIOPT.
Lasserre, J.B., Henrion, D., Prieur, C., and Trélat, E.
(2008). Nonlinear optimal control via occupation measures and lmi-relaxations. SICON.
Ledzewicz, U. and Schättler, H. (2007). Optimal controls for a model with pharmacokinetics maximizing
bone marrow in cancer chemotherapy. Mathematical
biosciences.
Lenhart, S. and Workman, J.T. (2007). Optimal control
applied to biological models. CRC Press.
Li, X., Omotere, O., Qian, L., and Dougherty, E.R. (2017).
Review of stochastic hybrid systems with applications
in biological systems modeling and analysis. EURASIP
Journal on Bioinformatics and Systems Biology.
Majumdar, A., Vasudevan, R., Tobenkin, M.M., and
Tedrake, R. (2014). Convex optimization of nonlinear
feedback controllers via occupation measures. IJRR.
Mevissen, M., Lasserre, J.B., and Henrion, D. (2011).
Moment and SDP relaxation techniques for smooth approximations of problems involving nonlinear differential equations. IFAC.
Mobilia, N. (2015).
Méthodologie semi-formelle
pour l’étude de systèmes biologiques: application à
l’homéostasie du fer. Ph.D. thesis, UGA.
Noel, V., Vakulenko, S., and Radulescu, O. (2011). Algorithm for identification of piecewise smooth hybrid
systems: application to eukaryotic cell cycle regulation.
In WABI.
Pace, A.M. and Burden, S.A. (2017).
Piecewisedifferentiable trajectory outcomes in mechanical systems subject to unilateral constraints. In HSCC.
Pakniyat, A. and Caines, P.E. (2014). On the relation
between the minimum principle and dynamic programming for hybrid systems. In CDC.
Rocca, A., Dang, T., Fanchon, E., and Moulis, J.M. (2016).
Application of the reachability analysis for the iron
homeostasis study. In HSB.
Rumschinski, P., Borchers, S., Bosio, S., Weismantel, R.,
and Findeisen, R. (2010). Set-base dynamical parameter estimation and model invalidation for biochemical
reaction networks. BMC systems biology.
Shia, V., Vasudevan, R., Bajcsy, R., and Tedrake, R.
(2014). Convex computation of the reachable set for
controlled polynomial hybrid systems. In CDC.
Soldatova, L.N., Aubrey, W., King, R.D., and Clare, A.
(2008). The exact description of biomedical protocols.
Bioinformatics.
Streif, S., Rumschinski, P., Henrion, D., and Findeisen,
R. (2013). Estimation of consistent parameter sets
for continuous-time nonlinear systems using occupation

measures and LMI relaxations. In CDC.
Vinter, R. (1993). Convex duality and nonlinear optimal
control. SICON.
Wardi, Y., Egerstedt, M., and Hale, M. (2015). Switchedmode systems: gradient-descent algorithms with armijo
step sizes. Discrete Event Dynamic Systems.
Xu, X. and Antsaklis, P.J. (2004). Optimal control
of switched systems based on parameterization of the
switching instants. IEEE trans. on automatic control.
Zhao, P., Mohan, S., and Vasudevan, R. (2017). Optimal
control for nonlinear hybrid systems via convex relaxations. arXiv preprint arXiv:1702.04310.
APPENDIX
The parameters of the haemoglobin production model are
shown in Table 5. These are taken from the column pmean
of (Bouchnita et al., 2016, Table 1). The ODEs used in the
hybrid system H are (fctrl ) and (frad ). The parameters
k2 and k7 are set to 0. The parameters ak3 and bk3 are
not used, as we search a new function to model k3 (t).
The initial condition of the controlled hybrid system, H,
is determined by the initial condition in Bouchnita et al.
(2016): (i0 , x(0)) = (1, (0, Fe0 , 0, 0, 0)).
Parameter
k1
k4
k5
k6
k8
Fe0
Feex
59Fe
ex

Numerical Value
10−3

1.410 ×
4.4710 × 10−4
7.2710 × 10−6
4.4710 × 10−4
1.1410 × 10−5
3.2110 × 102
4
3000

Units
s−1
s−1
fL.molecules−1 .s−1
s−1
s−1
atoms/fL
atoms/fL
atoms/fL

Table 5. The parameters (not scaled) for the
ODEs fctrl and frad , taken from the parameters
set pmean in Bouchnita et al. (2016).

x˙c




 x˙1
x˙2



x˙

 3
x˙4

= 1
= k1 Feex − u(t)x1 (t)
= − k4 x2 (t) − 4k5 x2 (t)x3 (t) + u(t)x1 (t)
= k6 x2 (t) − 4k5 x2 (t)x3 (t)
= k5 x2 (t)x3 (t) − k8 x4 (t)


x˙c




x
˙1




x
˙2





x
 ˙3
x˙4



x˙5




x˙6





x˙

 7
x˙8

=
=
=
=
=

1
k1 Feex − u(t)x1 (t)
− k4 x2 (t) − 4k5 x2 (t)x3 (t) + u(t)x1 (t)
k6 x2 (t) − 4k5 x2 (t)x3 (t)
k5 x2 (t)x3 (t) − k8 x4 (t)

= k1 59Feex − u(t)x5 (t)
= − k4 x6 (t) − 4k5 x6 (t)x7 (t) + u(t)x5 (t)
= k6 (x2 + x6 )(t) − 4k5 (x2 (t) + x6 )x7 (t)
= k5 x6 (t)x7 (t) − k8 x8 (t)

(fctrl )

(frad )

