Occupation measure methods for modelling and analysis of biological hybrid automata
Alexandre Rocca, Marcelo Forets, Victor Magron, Eric Fanchon, Thao Dang
OOccupation measure methods for modellingand analysis of biological hybrid systems
Alexandre Rocca ∗ , ∗∗ Marcelo Forets ∗ Victor Magron ∗ Eric Fanchon ∗∗ Thao Dang ∗∗ 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 Abstract:
Mechanistic models in biology often involve numerous parameters about which wedo not have direct experimental information. The traditional approach is to fit these parametersusing extensive numerical simulations (e.g. by the Monte-Carlo method), and eventually revisingthe model if the predictions do not correspond to the actual measurements. In this work wepropose a methodology for hybrid system model revision, when new types of functions are neededto capture time varying parameters. To this end, we formulate a hybrid optimal control problemwith intermediate points as successive infinite-dimensional linear programs (LP) on occupationmeasures. Then, these infinite-dimensional LPs are solved using a hierarchy of semidefiniterelaxations. The whole procedure is exposed on a recent model for haemoglobin production inerythrocytes.
Keywords: biological modelling, hybrid dynamical system, optimal control problem,semidefinite optimization, occupation measures.1. INTRODUCTION AND CONTEXTMechanistic models in biology generally involve many pa-rameters. The value of a given parameter can be eithermeasured directly in a dedicated experiment (e.g. mea-surement of a kinetic parameter of a biochemical reactionin enzymology), or inferred from data which provide rela-tionships between parameters and other biological entities.In this paper, we work with biological mechanisms mod-elled by ordinary differential equations (ODEs) and hybriddynamical systems. The motivation for using such modelsis to give quantitative predictions, when sufficient data isavailable to validate the model. And whenever new dataand knowledge of various types become available, they canbe incorporated within a formal framework.A basic issue in biological systems modelling is the deter-mination of numerical values for the parameters, or moregenerally a subset of the parameter space, under which themodel agrees to some extent with the available data. Wefocus on multiple-step experiments, in which a biologicalsystem is perturbed or measured during its evolution.In the biological systems modelling literature, it is com-mon to synthesise parameters using a Monte-Carlo sam-pling of the parameter space, which is validated then bynumerous simulations. An important effort to formalizeand validate the parameter synthesis of biological modelshas been made in works such as Donz´e (2010); Mobilia(2015); Dreossi (2016); Beneˇs et al. (2016), or Rumschinskiet al. (2010). Other such as Cardelli et al. (2017) or For all authors emails follow [email protected]
Bartocci et al. (2013) design ODE models satisfying setsof temporal constraints. When model simulation does notreproduce satisfactorily available experimental data, to adegree which depends on data quality, for any admissibleparameter value, the model has to be revised. One wayof revising the model is to represent parameters usingdifferent types of functions of time, reflecting underlyingbiological mechanisms. We introduce a systematic way,based on formal methods, to study mechanistic biologi-cal models in their experimental context and revise pa-rameters to produce conservative results with respect toexperimental data. In this work, we consider a problemof model revision, defined as finding time varying laws ofparameter evolution that minimizes the error in matchingexperimental measurements. Concretely, it is the followingoptimization problem:inf( x , u ) n exp (cid:88) j =1 dist( m ( x ( T j )) , z j ) (1)where x is a vector of biological variables, such as concen-trations, whose evolution is modelled by trajectories of abiological dynamical system. Time varying parameters arerepresented by the input variables u (modelling biologicalparameters) such that ∀ t ∈ [0 , T ] , u ( t ) ∈ U . X is the setof initial conditions and the set of pairs { ( T j , z j ) } j is theset of data points, for 1 ≤ j ≤ n exp , in the time frame[0 , T ]. An experimental measurement is a function of thevariables x and is modelled via the function m ( x ).The framework of our approach is a mathematical for-malization of experimental protocols as hybrid dynami-cal systems, describing biological systems of interest andexperiments which are performed on them. However, the a r X i v : . [ c s . S Y ] M a r lgorithm we provide can be applied to any biologicalhybrid system with similar model revision problems.In this paper we address the parameter synthesis andmodel revision problems (1) by formulating a particularinstance of the optimal control problem with intermediatecosts, when the objective function depends on the systemstrajectory and control inputs at a given set of timepoints. This problem is then approximated by multiplehybrid optimal control problems (HOCP) with one finalcost. To this end, we apply a recently developed methodof Zhao et al. (2017) from the field of certified convexoptimization to globally solve these HOCP. However, Themethod described in Zhao et al. (2017) produces piecewiseoptimal control functions which either may not correspondto biological knowledge of parameter variations or maybe difficult to yield coherent and meaningful biologicalinterpretations. Consequently, in order to respect realisticconstraints on parameters, we use smooth approximationsof the generated control input, in order to revise the givenmodel while maintaining good data fitting accuracy.The method is demonstrated on a hybrid system modellinghaemoglobin production presented in Section 4.1. The hybrid formalism has previously been used as an ab-straction method to simplify complex mechanisms whichare hard to analyse as seen in Noel et al. (2011); Roccaet al. (2016), or to represent “jump” evolution such as acti-vation processes in genes regulatory networks for exampleusing the stochastic formalism as in Li et al. (2017).Optimal control theory and variation theory have beenapplied to biological systems in several works. Most ofthem address the classical problem of finding a correctinput such that the system reaches a desired state. Forexample, one can control drug input such that a patientattains a healthy state, see Ledzewicz and Sch¨attler (2007)or Caraguel et al. (2016). Another example is the controlof some input in population studies as detailed in Bodineet al. (2008). A detailed review on the use of optimalcontrol in systems biology can be found in Lenhart andWorkman (2007). The problem of parameter estimationin 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 hy-brid 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); Xuand Antsaklis (2004); Bengea and DeCarlo (2005). Moregenerally, Pakniyat and Caines (2014) relies on DynamicProgramming and an extension of Pontryagin’s MaximumPrinciple. However, these approaches need a-priori knowl-edge either on the sequence of discrete transitions, or onthe number of visited subsystems. To perform optimalcontrol on hybrid systems, we build our work on the tech-niques from Zhao et al. (2017), which proposes a methodto obtain a global solution for hybrid systems with state-dependent transitions, without any a-priori knowledge onthe execution and the sequence of transitions. We refer to(Zhao et al., 2017, Section 1.1) and references therein formore details on optimal control of hybrid systems. Semidefinite programming (SDP) eases the resolution ofhard optimization problems and yields conservative re-sults ensured by positivity certificates. In Lasserre (2001),hierarchies of semidefinite relaxations were introducedfor static polynomial optimization. The definition of aninfinite-dimensional linear program (LP) over occupationmeasures, for optimal control problems, was first intro-duced in Vinter (1993). From this infinite-dimensional LP,Lasserre et al. (2008) defines hierarchies of Linear MatrixInequalities (LMI) relaxations, to synthesise a sequenceof polynomial controls converging to the solutions of theoptimal control problem. In Abdalmoaty et al. (2013) theauthors propose an extension to piecewise affine systems.Our underlying idea of constructing a suboptimal controlwith an iterative algorithm is similar to (Abdalmoatyet al., 2013, Section 4). However, we use this scheme tofind input functions allowing to reproduce data not onlyat a final time point but also at intermediate time points.We make use of the recent method proposed in Zhaoet al. (2017), which relies on occupation measures and asequence semidefinite relaxations to produce a sequence ofpolynomial controls converging to the optimal solution ofa HOCP. There exist other methods which use occupationmeasures and LMI relaxations to produce both admissiblecontrols and converging outer-approximations of the back-ward reachable set (BRS) Majumdar et al. (2014); Shiaet al. (2014), or the region of attraction (ROA) Kordaet al. (2014).Finally, we note that finding a sequence ofconverging outer-approximations for all valid parameterssets, such as in Streif et al. (2013); Dreossi (2016), isanother crucial issue in the context of systems biology.When dealing with hybrid systems, Shia et al. (2014) canbe applied to solve this problem.
Given a hybrid system modelling a multiple-step biolog-ical protocol and a set of experimental measurements,we propose a numerical algorithm to solve the modelrevision problem. This algorithm generates an admissibleinput function such that the revised model reproduces theexperimental measurements, and can be interpreted. Werelease a software package , in MATLAB, implementingthis algorithm. We evaluate this algorithm and its imple-mentation, on the haemoglobin production model takenfrom Bouchnita et al. (2016).The rest of the paper is organized as follows: in Section 2we give the necessary notations on dynamical hybrid sys-tems and the optimal control problem. Section 3 presentsour main contribution to the resolution of the hybridoptimal control problem with intermediate points. Finally,in Section 4 we study and discuss the model revision ofthe haemoglobin production model taken from Bouchnitaet al. (2016). 2. PRELIMINARIESWe first give the notations and basic notions on controlledhybrid systems, as well as the definition of the HOCP.Given x ∈ R n , let x i denote its i-th component. In general,letters in bold font denote multidimensional elements, and https://github.com/biosdp/biosdp ormal font unidimensional ones. Let B := { true , false } be the set of Booleans. Let R [ x ] denote the ring of realpolynomials in x ∈ R n , and let R d [ x ] be the subspace ofpolynomials whose degree is at most d . Let T be the timeinterval [0 , T ], where T is the final time (possibly ∞ ).Consider the n -dimensional ODE with inputs:˙ x ( t ) = f ( t, x ( t ) , u ( t )) , (2)with f : T × R n × R m → R n a vector field which is Lipschitzcontinuous in x and piecewise continuous in u . Let X and U be compact subsets of R n and R m respectively. Here, u : T → U is a feasible input function which representstime varying parameters, or external commands. The tuple F := ( T , X , U , f ) (3)defines a continuous dynamical system.We recall the definition of controlled hybrid systems adynamic hybrid systems formalism defined in Zhao et al.(2017). Definition 1. (Controlled hybrid system). A controlled hy-brid 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 n modes thenumber of modes. • E ⊆ I × I is the set of transitions e = ( i, j ) betweentwo modes: i is the source mode, and j the destinationmode. • X := (cid:96) i ∈I X i is the disjoint union of domains of H and X i the domain of the mode i . We note that X i is a compact subset of R n i with n i the dimensionof the mode i . The disjoint union (cid:96) can simplybe considered as a labelling operation on the set ofdomains by I , that is the set of mode indices. • U is the set of input values of H . • F := {F i } i ∈I is the set of continuous dynamicalsub-systems associated to each mode. The dynamicalsystem associated to mode i is: F i := ( T , X i , U , f i ) , with f i : T × X i × U → R n i a vector field polynomialin x and affine in u . • S := (cid:96) e ∈E S e is the disjoint union of guards S e ⊆ X i associated to each transition e = ( i, j ) ∈ E . The guard S ( i,j ) defines the switch condition from i to j : for x ∈ X i , if x ∈ S ( i,j ) then the system at x can makethe transition from mode i to mode j . • R := {R e } e ∈E is the set of reset maps, each resetmap R e : S e → X j being associated to a transition e := ( i, j ) ∈ E and it defines how the continuousvariables may change after the discrete transitionfrom mode i to mode j .Additionally, the CHS defined in Zhao et al. (2017) mustrespect a few assumptions. All the guards S ( i, · ) are disjoint,and S ( i,j ) ⊆ ∂ X i , with ∂ X i designating the border of X i ,for each pair of modes i and j . The initial set is restrictedto a single point x , with an associated mode i . Thevector fields f i are affine in u and have a nonzero normalcomponent on the boundary of X i . These assumptionsensure that any CHS is deterministic. Noting λ ( x ( t )) thefunction which associates to an instantaneous state x ( t )its corresponding mode, these assumptions ensure that themode corresponding to x ( t ) is unique. Given a CHS H , the hybrid optimal control problem isdefined as follows. Let X , and X T , be the initial set, andtarget set defined by: X := (cid:97) i ∈I X ,i , X T := (cid:97) i ∈I X T,i , (4)where X i and X T i are a compact subsets of X i foreach mode i ∈ I . Let i and i T be the initial modeand the final mode at time T , respectively. Then, given( i , x (0)) = ( i , x ) ∈ X and u : T → U an inputfunction, we say that for T >
0, ( x ( t ) , u ( t )) ∈ P is anadmissible pair on T and P is the set of admissible pairs,if ( i, x ( t )) ∈ X is a trajectory of CHS as defined by (Zhaoet al., 2017, Algorithm 1) and ( i T , x ( T )) ∈ X T . Finally,the hybrid optimal control problem for a CHS H , (HOCP),is defined by: J ∗ hocp := inf ( x , u ) ∈P (cid:90) T h λ ( x ( t )) (cid:0) t, x λ ( x ( t )) ( t ) , u ( t ) (cid:1) dt + H λ ( x ( T )) (cid:0) x λ ( x ( T )) ( T ) (cid:1) . (5)where { h i : [0 , T ] × R n i × R m → R } i ∈I and { H i : R n i → R } i ∈I are set of measurable functions.3. OPTIMAL CONTROL FOR MODEL REVISIONIn our biological context hybrid dynamical systems canmodel multiple steps experiments. Indeed, experimentalprotocols associated to these experiment can be consideredas a set of concurrent processes where each process ismodelled as a hybrid system. Thus, the hybrid system ofthe protocol can be described as the parallel compositionof the hybrid system describing each process. It is thencrucial to proceed to model revision while taking intoaccount the biological system in the evolving environmentof the complete protocol. For this purpose, we propose inthis section a method to produce time varying parametersreproducing multiple data in the context of a multiple-phase protocol modelled by a hybrid system. However, wealso argue that this method can be use for more generalbiological systems modelled as CHS. We solve the modelrevision problem of a dynamical hybrid system modellinga biological system together with a set of experiments.Therefore, we search for parameters as time varying func-tions fitting a set of data points, defined in the introductionas in (1). In this aim, Section 3.1 first formulates (1) asa particular instance of the optimal control problem onhybrid systems with intermediate costs. Then, we proposea first approximation as a set of instances of the HOCP(5) defined previously. The solution of each sub problem isobtained using the previous results from (Zhao et al., 2017,Section 4). Finally, in Section 3.2 we explain the completealgorithm addressing our initial problem. We provide a method to find time varying parametersof biological system, modelled as input functions u ( t ), inorder to fit the hybrid systems model to a set of experi-mental data. Thus, we write our problem as an optimalcontrol problem where desired input functions are theoptimal controls which minimize the distance of the resultsproduced by the model and these experimental data. Ex-perimental measurements, represented by a function m ( x ),re performed at given specific times T j , 1 ≤ j ≤ n exp .Let z j be the observed value of the experimental measure-ment at time T j , then n exp is the number of experimentaldata points. Let X T j,i be compact subsets of X i , and X T j := (cid:96) i ∈I X T j ,i . As in (5), let ( i , x (0)) ∈ X , andsuppose that we are given a set of time values { T j } , with1 ≤ j ≤ n exp , and T n exp = T . Given an input function u : T → U , we say that ( x , u ) ∈ P int is an admissiblepair for a problem with intermediate points and P int theassociated set of admissible pairs, if ( i ( t ) , x ( t )) ∈ X is aCHS trajectory, and ( i T j , x ( T j )) ∈ X T j for all j .Let H ( x ( T j )) be a cost at time T j , and h ( t, x ( t ) , u ( t ))a running cost for the whole [0 , T ] interval. The optimalcontrol problem with intermediate points for the CHS H is then: J ∗ hocp := inf ( x , u ) ∈P int (cid:90) T h ( t, x ( t ) , u ( t )) dt + (cid:88) ≤ j ≤ n exp H ( x ( T j )) (6)In our biological context, H ( x ( T j )) = || m ( x ( T j )) − z j || . To our best knowledge there is no method to efficiently ad-dress directly problem (6), and obtain converging sequenceof solution in a similar manner to the simpler problem(5) address in Zhao et al. (2017). Thus, we search foran admissible solution using a greedy algorithm. Conse-quently, this solution is a good trade-off between compu-tation cost, optimality and flexibility as we can handlea large class of models. Moreover, this method does notconstraint the form of the sought times parameters: thisreduces the assumptions during the modelling and easesthe final biological interpretations, unlike methods basedon simulations. Given 1 ≤ j ≤ n exp , let J j ( t, x ( t ) , u ( t )) := (cid:90) T j T j − h ( t, x ( t ) , u ( t )) + H ( x ( T j )) , with T = 0, and T n exp = T , such that J ( t, x ( t ) , u ( t )) = (cid:88) ≤ j ≤ n exp J j ( t, x ( t ) , u ( t )) . Noting ( i ( j ) ( t ) , x ( j ) ( t )) a trajectory of a CHS H on theinterval T j := [ T j − , T j ], and similarly ˜ u ( j ) ( t ) the controlon T j , we consider the following problem as particularinstance of (5): J ∗ j := inf ( x ( j ) , ˜ u ( j ) ) J j ( t, x ( j ) ( t ) , ˜ u ( j ) ( t ))s.t. x ( j ) cont. part of a trajectory ( i ( j ) , x ( j ) ) on T j , ˜ u ( j ) ( t ) ∈ U , ∀ t ∈ T j , ( i ( j ) ( t ) , x ( j ) ( t )) ∈ X , ∀ t ∈ T j , if j = 1 , ( i (1) (0) , x (1) (0)) ∈ X , if j ≥ . ( i ( j ) ( T j − ) , x ( j ) ( T j − )) = ( i ( j − ( T j − ) , x ( j − ( T j − )) , ( i ( j ) ( T j ) , x ( j ) ( T j )) ∈ X T j . (7)We note that if a transition i → i (cid:48) occurs at the time T j ofthe interval [ T j − , T j ], we retain only the left part in the mode i for the next optimization on the interval [ T j , T j +1 ].Let ˜ u ( t ) and ( i ( t ) , x ( t )) be respectively the control andthe trajectory, for t ∈ [0 , T ]. They are respectively definedby the concatenation of all the controls ˜ u ( j ) ( t ) and thetrajectories ( i ( j ) ( t ) , x ( j ) ( t )) on the sub-intervals [ T j − , T j ].By construction, ( x ( t ) , ˜ u ( t )) is an admissible pair for (6),as ( i T j , x ( T j )) = ( i ( j ) T j , x ( j ) ( T j )) ∈ X T j . Remark 1.
We emphasize that ( x ( t ) , ˜ u ( t )) is not necessaryan optimal solution for (6). Moreover, as the optimizationproblem (7) is obtained through a greedy scheme, we haveno guarantee that its optimal cost J ∗ j is inferior to a given ε . However, as we search to equally fit all the data pointssearching iteratively for the control is satisfiable solution. Let ( T j , z j ), 0 ≤ j ≤ n exp be pairs of experimental datapoints and their measurement time, and we also note i , and x the initial mode and initial conditions of thestudied CHS H respectively. Let r be a given startingrelaxation degree. Algorithm 1 finds an admissible solutionto (6), by solving a sequence of the HOCP (7) for eachexperimental data point ( T j , z j ). For each j , the degree ofthe polynomial control ˜ u ( j ) ( x ( t ) , t ) is determined as thesmallest degree such that || m ( x ( T j )) − z j || ≤ ε . Indeed,in the context of biological system modelling we desire toobtain a control of degree as small as possible to avoidoverfitting. We note that as explained in Remark 1, wecannot ensure the converge to a solution of accuracy ε .Thus, ε is only a stopping criteria. Then, for each iterationover j , Algorithm 1 is decomposed in three steps.The first step is the procedure HOCP , associated to aninstance of the HOCP (5) for j -th pairs ( T j , z j ). Givena relaxation order d r ≥ r , we solve the relaxed primaldefined in (Zhao et al., 2017, Section 5.1). From Zhao et al.(2017) method we first obtain M d r ( y µ i ), the sequence mo-ment matrices of degree d r , associated to the occupationmeasure µ i of each mode i ∈ I . We also obtain J ( d r ) j anunder approximation of the optimum of (7). The secondstep is the procedure Synth , which returns the admissiblecontrol ˜ u ( j ) ( t, x ) of degree d u ≤ d r using a truncatedmoment matrix M d u ( y µ i ) of M d r ( y µ i ) at the reduceddegree d u . The third and last step is the procedure Simu .It performs the validation that the synthesised control ˜ u ( j ) yields || m ( x ( T j )) − z j || ≤ ε . This step is done by approxi-mating the trajectory of the controlled hybrid system usinga solver of ODE with discrete events to produce numericalsimulations. If in iteration j , || m ( x ( T j )) − z j || ≤ ε , then x ( j ) ( T j ) and the corresponding mode i f reached at t = T j by the numerical simulations are the initial conditions forthe next iteration j + 1. Otherwise, the Ctrl Synth and
Simulate procedures are repeated while increasing thedegree of the synthesised polynomial control until d u = d r .In case the condition || m ( x ( T j )) − z j || ≤ ε is still notsatisfied, the relaxation order d r is increased, and the threesteps are repeated.If ε ≤ J ( d r ) j then we are sure that we the given initial con-dition at step j , there is no control such that || m ( x ( T j )) − z j || ≤ ε . Consequently, we keep our previous result ˜ u ( j ) and the corresponding mode i f reached at t = T j by theumerical simulations are the initial conditions for thenext iteration j + 1. Algorithm 1. procedure Algorithm H , { ( T j , z j ) } j , i , x , ε, r ) T init = 0 for all experimental data ( T j , z j ) do d u = , d r = r, err = + ∞ while err ≥ ε ∧ J ( d r ) j ≤ ε do J ( d r ) j , M d r ( y µ ) = HOCP ( H , i , ... ... x , T init , T j , z j , d r ) while err ≥ ε and d u ≤ d r do ˜ u ( j ) ( x ( t ) , t ) = Synth ( M d r ( y µ ) , d u ) ( i f , x ( j ) ( t )) = Simu ( H , ˜ u ( j ) ( x ( t ) , t ) , ... ... i , x , T init , T j ) err = H ( x ( j ) ( T j ) , z j ) increase d u end while increase d r end while i = i f x = x ( j ) ( T j ) T init = T j end for end procedure
4. RESULTS AND DISCUSSIONIn this section, using the method developed in Section 3,we revise the model of haemoglobin production by findinga better fit for the time varying parameter noted k inBouchnita et al. (2016), with respect to the same errorfunction. To be consistent with our notations from Sections2 and 3, we will note u ( t ) the input modelling the timevarying parameter k . The paper Bouchnita et al. (2016) addressed, on an exper-imental protocol model, the problem of refining the pa-rameters space, and to fit multiple data sets of experimen-tal results. In Bouchnita et al. (2016), the experimentalprotocol was not explicitly formalized as a hybrid system.However, we use the parameters value from this previouswork as starting instantiation to provide less restrictiveconstraints on the time varying parameter k .The ODEs ( f ctrl ) model the evolution of the haemoglobinproduction in the differentiating erythrocyte cells situatedin the bone marrow. The variables x to x representrespectively the internal iron in the cell Fe, the hemeH , the globin G, and the haemoglobin Hb. The hybridsystem H models an experimental protocol designed tomeasure the integration of iron inside heme (H) at multiplemilestones of the cell differentiation. For example the datapoint at time t = 7 hours in Table 3, is obtained throughthe following procedure: we first start with a control batchof cells, then at time t = 4 hours after the start ofthe experiment, the culture medium is perturbed with aninjection of measurable radioactive iron Fe for a subsetof the cells. This perturbation implies new ODEs ( f rad ) Table 1. Dimensions (with x c ), vector fields,domains, and input sets for the controlled hy-brid system H of the haemoglobin productionmodel. Mode n i f i ( t, x , u ) X i U i i = 1 5 f ctrl ( t, x , u ) [0 , × [0 , [0 , i = 2 9 f rad ( t, x , u ) [4 , × [0 , [0 , i = 3 5 f ctrl ( t, x , u ) [7 , × [0 , [0 , i = ... ... ... ... ... i = 13 5 f ctrl ( t, x , u ) [45 , × [0 , [0 , i = 14 9 f rad ( t, x , u ) [52 , × [0 , [0 , e = ( i, j ) S e R e i = 1 (1 , t == 4 (cid:104) I , O , (cid:105) i = 2 (2 , t == 7 (cid:2) I , , O , (cid:3) i = 3 (3 , t == 8 (cid:104) I , O , (cid:105) i = ... ... ... ... i = 13 (13 , t == 52 (cid:104) I , O , (cid:105) Table 2. Transitions, guards, and reset maps ofthe controlled hybrid system H .modelling the evolution of two inter-dependant models: themodel of non-radioactive haemoglobin production and themodel of haemoglobin production with radioactive species.Three hours after the perturbation with radioactive iron,the total radioactive heme is measured, meaning the hemefree in the cell and the one in the radioactive haemoglobin.This measurement is given by the formulas H + 4 Hb.The data in Table 3 is then the observed radioactivity di-vided by three hours. Finally, these measurements provideresults on the variation during the cell differentiation ofthe integration of iron in heme, which is associated to theparameter k .The controlled hybrid system H associated to the ex-periment studied in Bouchnita et al. (2016) and thehaemoglobin production model are given, in a shortenedversion, in Tables 1 and 2. The ODEs ( f ctrl ) and ( f rad ), aswell as the numerical values of the parameters, are givenin Appendix. In the implementation, we also introduce avariable x c modelling time, whose derivative is equal to 1.For numerical reasons, it is necessary to scale the param-eters and state variables, making it easier for the solverto succeed in solving the relaxed problem. Similarly, tofacilitate the numerical optimization we rewrite the controlvariable u ( t ) ∈ U = [0 ,
1] as u ( t ) = ζ ˆ u ( t ), with ζ (cid:28) u ( t ) ∈ [0 , /ζ ]. While the scale factor ζ may take differentvalues depending on the numerical optimization details,the objective control u ( t ) always evolves in [0 , H , wesolve the optimal control problem with intermediate timepoints defined in (6), using the method from Section 3.1 Time [hours] k ( t ) [ hou r s - ] × -4 Generated optimal control and associated approximating functions
Generated Optimal control u gen (t) Piecewise polynomial u poly (t) Hill function u hill (t) Step function u step (t)
Fig. 1. Synthesised optimal control and various approxi-mations that yield a realistic interpretation.
Time [hours] F e59 ( t ) [ m o l / m ] Fe59(t)
Time [hours] G t o t ( t ) [ m o l / m ] Gtot(t)
Time [hours] H ( t ) + * H b59 ( t ) [ c p m / ( - L ) h - ] Observed variable: H59(t)+4*Hb59(t)
Fig. 2. Radio-active variables Fe, G tot corresponding to x and x in ( f rad ), as well as, the comparison of themeasurement function results to the data.and its implementation in Section 3.2. The experimentalmeasurement is modelled by the function m ( x ) = x +4 x .Thus, we set H ( x ( T j )) := ( x ( T j ) + 4 x ( T j ) − z j ) , as wesearch to minimize the total residual error term: ε total = (cid:88) ≤ j ≤ n exp (cid:112) H ( x ( T j )) (cid:80) ≤ j ≤ n exp z j . (8)The original experimental data points ( T j , z j ) are given inTable 3. Here, the input control u ( t ) models some hidden Time ( h ) 7 11 19 27 35 45 55Measure ( cpm e − L · h − ) 16 85 348 391 399 481 395 Table 3. Experimental data points ( T j , z j ) usedas references.mechanism which evolves with the differentiation of thecells. It should be the same function of time for both thecontrol and the radioactive cells batch. However, as thecontrol generated by Algorithm 1 is piecewise for eachmode, and the fact that our data are on the radioactivespecies only, the solution of the optimization problem withonly a final cost H ( x ( T i )) is not balanced , having a muchstronger control in the modes where the radioactive speciesare evolving.A workaround for the balancing problem is the following.We add a small penalization cost c i ( t ) = (0 . u ( t )) toequilibrate the control when i corresponds to a modewith radioactive species, otherwise c i ( t ) = 0. In a similar vein, we add another penalization cost c i ( t ) = ( u ( T j ) − u ( t )) to avoid when the control strongly varies betweentwo iterations j on the interval [ T j − , T j ] and j + 1 on[ T j , T j +1 ] (with the exception of the first iteration). Thisleads to h i ( t, x ( t ) , u ( t )) = c i ( t ) + c i ( t ). Let us note that,even if these additional costs can eventually degrade theaccuracy of the data fitting, we gain in terms of biologicalinterpretation of the resulting traces.Finally, by partitioning the computation in the time do-main, we can greatly reduce the computational cost ateach iteration. More technically, since the transitions ofthe hybrid system H are fully determined by the time t , we can pre-compute the function λ : R + → I , whichassociates a mode λ ( t ) to each time instant t . Thus, eachiteration j of Algorithm 1 can be restrained to the sub-hybrid system H j of H , constituted by the modes visitedin the interval [ T j − , T j ].For numerical implementation, the problem on measuresis formulated in SPOTLESS , and then we extract theprimal solution provided by a primal-dual SDP solver. Todo so, we use the implementation from Zhao et al. (2017)to generate the dual problem. As an SDP solver we usedMOSEK Andersen and Andersen (2000) v.7.1. These toolsare used in MATLAB v.9.0 (R2016a). Performance resultsare obtained with an Intel Core i7-5600U CPU (2 .
60 GHz)with 16Gb of RAM running on Debian 8. We only solvethe problem for a relaxation order r = 4, as any higherorder would be too memory expensive. In this particularexample as constant input leaded to satisfying results, wedid not impose any constraint (cid:15) in the algorithm. Usingthis configuration, the total time taken by Algorithm 1is 2107s, with 1700s spent in the HOCP procedure, and390s in the
Synth procedure. On Figure 1, the controlgenerated by Algorithm 1 is shown in blue. This controlis piecewise, and clearly divided in two phases: before andafter t equals 11 hours. However, the control synthesisedis still difficult to interpret as a biological phenomenon.Consequently, we propose three additional fits of thiscontrol to ease interpretation by using functions closer tobiological knowledge. In Table 4 one can find the totalerror associated to all the possible controls, as well as theprevious result of Bouchnita et al. (2016). In Figure 2,we show a graphical representation of how closely eachfunction can control the model to reach the desired datapoints. Control Type ε total Original Control Bouchnita et al. (2016) 0 . . . . . Table 4. Total error ε total associated to eachproposed input. The SPOTLESS implementation was taken fromhttps://github.com/spot-toolbox/spotless .2 Discussion
In a simulation-based approach, we have to propose forthe desired time varying parameter, a template functionto fit the data, e.g. a polynomial of given degree. If wewant to fit a polynomial of higher degree, the simulationshave to be run again multiple times. On the contrary, theproposed approach returns a control signal, and since thefit to data points is performed a posteriori, there is noadditional computation cost in refining the model.From the form of the experimental data points, an usualhypothesis is that u ( t ) should be similar to a jump func-tion, with a low value for the two first points, and ahigher one for the following ones. However, even withsuch information a good fit is not easily achieved withsimulations.The control generated with Algorithm 1 returns the ex-pected “jump” behaviour for u ( t ), and even with a lowrelaxation degree, the total residual error for the generatedcontrol is 9 .
59% which is much lower than the 22 .
8% fromBouchnita et al. (2016).We first fit a step function to the generated control, witha change at t = 11. The associated error of 12 .
24% isstill lower than Bouchnita et al. (2016), yet being higherthan the generated control mainly due to the second-to-last point.The second fit is a piecewise polynomial function in twopieces. The first piece, for t ∈ [0 , t ∈ [11 , . a priori on a particular form, avoiding the need forextensive numerical simulations. The generated control isaccurate, and computed in a reasonable time ( ∼ ECC .Andersen, E.D. and Andersen, K.D. (2000). The MosekInterior Point Optimizer for Linear Programming: AnImplementation of the Homogeneous Algorithm.Bartocci, E., Bortolussi, L., and Nenzi, L. (2013). Atemporal logic approach to modular design of syntheticbiological circuits. In
CMSB .Beneˇs, N., Brim, L., Demko, M., Pastva, S., and ˇSafr´anek,D. (2016). Parallel smt-based parameter synthesis withapplication to piecewise multi-affine systems. In
ATVA .Bengea, S.C. and DeCarlo, R.A. (2005). Optimal controlof switching systems. automatica .Bodine, E.N., Gross, L.J., and Lenhart, S. (2008). Optimalcontrol 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 oferythropoiesis and hemoglobin production.
JIOPM .Caraguel, F., Lesart, A.C., Est`eve, F., van der Sanden,B., and St´ephanou, A. (2016). Towards the design ofa patient-specific virtual tumour.
Computational andmathematical methods in medicine .ardelli, L., ˇCeˇska, M., Fr¨anzle, M., Kwiatkowska, M.,Laurenti, L., Paoletti, N., and Whitby, M. (2017).Syntax-guided optimal synthesis for chemical reactionnetworks. In
CAV .Donz´e, A. (2010). Breach, a toolbox for verification andparameter synthesis of hybrid systems. In
CAV .Dreossi, T. (2016). Sapo: Reachability computation andparameter synthesis of polynomial dynamical systems. arXiv .Korda, M., Henrion, D., and Jones, C.N. (2014). Controllerdesign and region of attraction estimation for nonlineardynamical systems.
IFAC .Lasserre, J.B. (2001). Global optimization with polynomi-als and the problem of moments.
SIOPT .Lasserre, J.B., Henrion, D., Prieur, C., and Tr´elat, E.(2008). Nonlinear optimal control via occupation mea-sures and lmi-relaxations.
SICON .Ledzewicz, U. and Sch¨attler, H. (2007). Optimal con-trols for a model with pharmacokinetics maximizingbone marrow in cancer chemotherapy.
Mathematicalbiosciences .Lenhart, S. and Workman, J.T. (2007).
Optimal controlapplied to biological models . CRC Press.Li, X., Omotere, O., Qian, L., and Dougherty, E.R. (2017).Review of stochastic hybrid systems with applicationsin biological systems modeling and analysis.
EURASIPJournal on Bioinformatics and Systems Biology .Majumdar, A., Vasudevan, R., Tobenkin, M.M., andTedrake, R. (2014). Convex optimization of nonlinearfeedback controllers via occupation measures.
IJRR .Mevissen, M., Lasserre, J.B., and Henrion, D. (2011).Moment and SDP relaxation techniques for smooth ap-proximations of problems involving nonlinear differen-tial equations.
IFAC .Mobilia, N. (2015).
M´ethodologie semi-formellepour l’´etude de syst`emes biologiques: application `al’hom´eostasie du fer . Ph.D. thesis, UGA.Noel, V., Vakulenko, S., and Radulescu, O. (2011). Al-gorithm for identification of piecewise smooth hybridsystems: application to eukaryotic cell cycle regulation.In
WABI .Pace, A.M. and Burden, S.A. (2017). Piecewise-differentiable trajectory outcomes in mechanical sys-tems subject to unilateral constraints. In
HSCC .Pakniyat, A. and Caines, P.E. (2014). On the relationbetween the minimum principle and dynamic program-ming for hybrid systems. In
CDC .Rocca, A., Dang, T., Fanchon, E., and Moulis, J.M. (2016).Application of the reachability analysis for the ironhomeostasis study. In
HSB .Rumschinski, P., Borchers, S., Bosio, S., Weismantel, R.,and Findeisen, R. (2010). Set-base dynamical param-eter estimation and model invalidation for biochemicalreaction networks.
BMC systems biology .Shia, V., Vasudevan, R., Bajcsy, R., and Tedrake, R.(2014). Convex computation of the reachable set forcontrolled 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 setsfor continuous-time nonlinear systems using occupation measures and LMI relaxations. In
CDC .Vinter, R. (1993). Convex duality and nonlinear optimalcontrol.
SICON .Wardi, Y., Egerstedt, M., and Hale, M. (2015). Switched-mode systems: gradient-descent algorithms with armijostep sizes.
Discrete Event Dynamic Systems .Xu, X. and Antsaklis, P.J. (2004). Optimal controlof switched systems based on parameterization of theswitching instants.
IEEE trans. on automatic control .Zhao, P., Mohan, S., and Vasudevan, R. (2017). Optimalcontrol for nonlinear hybrid systems via convex relax-ations. arXiv preprint arXiv:1702.04310 .APPENDIXThe parameters of the haemoglobin production model areshown in Table 5. These are taken from the column p mean of (Bouchnita et al., 2016, Table 1). The ODEs used in thehybrid system H are ( f ctrl ) and ( f rad ). The parameters k and k are set to 0. The parameters a k and b k arenot used, as we search a new function to model k ( t ).The initial condition of the controlled hybrid system, H ,is determined by the initial condition in Bouchnita et al.(2016): ( i , x (0)) = (1 , (0 , Fe , , , Parameter Numerical Value Units k . × − s − k . × − s − k . × − fL.molecules − . s − k . × − s − k . × − s − Fe . × atoms/fLFe ex Fe ex Table 5. The parameters (not scaled) for theODEs f ctrl and f rad , taken from the parametersset pmean in Bouchnita et al. (2016). ˙ x c = 1˙ x = k Fe ex − u ( t ) x ( t )˙ x = − k x ( t ) − k x ( t ) x ( t ) + u ( t ) x ( t )˙ x = k x ( t ) − k x ( t ) x ( t )˙ x = k x ( t ) x ( t ) − k x ( t ) ( f ctrl ) ˙ x c = 1˙ x = k Fe ex − u ( t ) x ( t )˙ x = − k x ( t ) − k x ( t ) x ( t ) + u ( t ) x ( t )˙ x = k x ( t ) − k x ( t ) x ( t )˙ x = k x ( t ) x ( t ) − k x ( t )˙ x = k Fe ex − u ( t ) x ( t )˙ x = − k x ( t ) − k x ( t ) x ( t ) + u ( t ) x ( t )˙ x = k ( x + x )( t ) − k ( x ( t ) + x ) x ( t )˙ x = k x ( t ) x ( t ) − k x ( t ) ( f radrad