Parameter inference in mechanistic models of coupled differential equations is a challenging problem. We propose a new method using kernel ridge regression in Reproducing Kernel Hilbert Spaces (RKHS). A three-step gradient matching algorithm is developed and applied to a realistic biochemical model

Filippone, Maurizio

Husmeier, Dirk

Niu, Mu

Rogers, Simon

English

Enlighten

Inference in nonlinear differential equationsMu Niu1, Maurizio Filippone2, Dirk Husmeier1, Simon Rogers21 School of Mathematics and Statistics, University of Glasgow, Glasgow UK2 School of Computing Science, University of Glasgow, Glasgow UKE-mail for correspondence: mu.niu@glasgow.ac.ukAbstract: Parameter inference in mechanistic models of coupled differentialequations is a challenging problem. We propose a new method using kernel ridgeregression in Reproducing Kernel Hilbert Spaces (RKHS). A three-step gradientmatching algorithm is developed and applied to a realistic biochemical model.Keywords: RKHS; nonlinear ordinary differential equation; gradient matching.1 IntroductionMany processes in science and engineering can be described by dynamicalsystems models based on nonlinear ordinary differential equations (NODEs).However, the parameters of the NODEs are often unknown and not directlymeasurable. Direct inference requires a computationally expensive numer-ical integration of the NODEs every time the parameters are adapted.For that reason, approximate methods based on gradient matching haverecently gained much attention; see e.g. Dondelinger et al. (2013) andHeinonen et al.(2014). The purpose of the present article is to try a newvariant of this approach based on reproducing kernel Hilbert spaces (RKHS).We consider systems governed by first order multivariate NODEs:x˙ =dxdt= f(x(t), θ), (1)with initial value x(t1) = x1, where x(t) is an r dimensional vector of statevariables (x(t) ∈ Rr). We observe the states of the system at n time points(y1, ..., yn) and assume that the observations consist of the states corruptedby Normal additive iid noise yi = x(ti) + i, where i ∼ N(0, σ2) iid. Inthis work, a framework of penalized regression for vector-valued functionsin RKHS is used to make predictions of the system state. A three stepgradient matching approach is developed to learn the NODEs parameters:This paper was published as a part of the proceedings of the 30th Interna-tional Workshop on Statistical Modelling, Johannes Kepler Universita¨t Linz, 6–10July 2015. The copyright remains with the author(s). Permission to reproduce orextract any parts of this abstract should be requested from the author(s).2 Inference in nonlinear differential equations1. Learn a smooth function g over the observations using an RKHSapproximation with observation (y1, ..., yn) at time points (t1, ..., tn).2. Given g, learn the NODEs parameters θ by gradient matching tominimize the difference between g˙(t) and f(g(t), θ).3. Re-estimate g using the optimised θ and observed trajectory y. Givenre-estimated g, optimise θ as in step 2.2 MethodologyFor each state variable xs indexed by s = 1, ..., r, a positive definite kernelks : R × R → R is used to define the Hilbert space Gs of the smootherfunction gs. The squared exponential kernel (Gaussian kernel) is used inthis study with lengthscale parameter l: k(tk, ti) = exp( − l−2(tk − ti)2).In step one, kernel ridge regression is used to learn gs, the smoother of thesth state with the following loss function:L1(ls) =n∑i=1(gs(ti)− ysi)2 + ||gs||2 (2)gs(ti) =n∑k=1bskks(tk, ti) (3)bs = (Ks + λsI)−1ys, ||gs||2 = λsbTsKsbs, (4)where ks(., ti) is the ith basis function and Ks is the Gram kernel matrix. ysis the vector of observations for the sth state, ||gs|| is the norm of Gs and theregularisation factor λs is estimated using leave-one-out cross validation. Iis the n× n identity matrix. Different states are assumed to have differentlengthscales, ls, and these are estimated independently through a gradient-based quasi-Newton optimisation routine. With the optimised ls, the stateestimation gs(tj) of xs(tj) is made on a uniform grid of m time pointsindexed by j = 1, ...,m. The gradient of the smoother is:g˙s(tj) =n∑i=1bsidk(ti, tj)dti. (5)In step 2, the NODE parameters are estimated by gradient matching forall states with the loss function:L2(θ) =r∑s=1m∑j=1(g˙s(tj)− f(gs(tj), θ))2, (6)where f(gs(tj), θ)is the target gradient generated by feeding the outputof the sth smoother, gs(t), into the NODEs system. θ is optimised using agradient-based quasi-Newton routine.Finally, in step 3, we use the estimated θ from step 2 to construct the NODEsystem. The lengthscale vector l for all state variables is then re-estimatedNiu et al. 3using the full loss function.L3(l) =r∑s=1( n∑i=1(gs(ti)− ysi)2 + ||gs||2 +m∑j=1(g˙s(tj)− f(gs(tj), θ))2)(7)The updated smoother g˜(t) is calculated using l. The refined θ can then beoptimised iteratively by repeating the second step.3 Application: Calcium modelThe calcium model (Peifer and Timmer, 2007) represents the oscillationsof calcium signaling in eukaryotic cells via a dynamic system with statescorresponding to the concentrations of free calcium in cytoplasm Cac andendoplasmic reticulum Car, as well as active Gα. and phospholipase-C, Pc.dGαdt= k1 + k2Gα − k3PcR1(Gα)− k4CacR2(Gα)dPcdt= k5Gα − k6R3(Pc)dCacdt= k7PcCacR4(Car) + k8Pc + k9Gα − k10R5(Cac)− k11R6(Cac)dCardt= −k7PcCacR4(Car) + k11R6(Cac) (8)where Ri(x) =xx+kmi. We followed Oates et al.(2014) and fixed the kmiparameters, leaving k1:11 to be inferred. Three scenarios were tested: thenoise-free case, and adding iid Gaussian noise with signal-to-noise ratioSNR = 50db and SNR = 10db to the numerical solution of the model.We sampled n = 100 (t1 = 0, t100 = 20) regularly spaced observations, andmade state predictions on m = 200 grid points. For 10db and 50db noise,100 independent data sets were generated.The mean NODEs parameter estimates with error bars are shown in Fig-ure 1. We see that when no noise is added, our estimates closely agree withthe true values. The parameter estimates at SNR = 50db noise are similarto the noise-free case when the parameters are smaller than 10, but driftaway from the true value for parameters larger than 10. For the highestnoise level, SNR = 10db, the NODEs parameter estimates become worse,although we still have 4 reliable estimates out of 11 parameters. For theSNR = 50db case, the length scale parameters l1:4 for the 4 states at theoptimum after step 2 are [0.06, 0.18, 0.01, 0.05]. After step 3, these change to[0.06,0.81,0.22,0.79], suggesting that a regularising effect from the NODEsexists.4 ConclusionWe have described an RKHS based gradient matching approach for pa-rameter inference in a system of NODEs. In low noise scenarios we obtain4 Inference in nonlinear differential equations0 1 2 3 4 5 6 701234567llllllllno noise50db10db0 50 100 150050100150llllno noise50db10db10 25102550100150FIGURE 1: Scatter plot of true versus estimated NODEs parameters in threenoise scenarios. No noise (black crosses), SNR = 50db (red circles) and SNR =10db (blue triangles) noise. Left and right plots show parameters with true values< 10 and > 10 respectively.good parameter estimates for a realistic biochemical application, and wehave quantified the deterioration resulting from increased noise levels. Wehave demonstrated that in the proposed three-step algorithm, the lengthscale parameters of the kernels are regularised by the gradient matchingstep. A potential explanation for the poor performance in the high-noiseregime is related to the choice of kernel. In our future work, we will ex-plore the effect of using non-stationary kernels (like the MLP kernel) andless smooth kernels (like the Mate´rn class kernels) as alternatives to thesquared exponential kernel.Acknowledgments: This project is funded by the EPSRC (EP/L020319/1).ReferencesDondelinger, F. et al. (2013) ODE parameter inference using adaptive gra-dient matching with Gaussian processes. JMLR W&CP 31: 216-228.Heinonen, M., and d’Alch-Buc, F. (2014). Learning nonparametric differ-ential equations with operator-valued kernels and gradient matching.arXiv preprint arXiv:1411.5172.Oates, C.J., et al. (2014). Causal network inference using biochemical ki-netics. Bioinformatics,30(17), i468-i474.Peifer, M. and Timmer, J. (2007). Parameter estimation in ordinary dif-ferential equations for biochemical processes using the method ofmultiple shooting. Systems Biology, IET,1.2., 78 – 88.

Inference in Nonlinear Differential Equations

http://eprints.gla.ac.uk/108342/1/108342.pdf

Inference in Nonlinear Differential Equations

Abstract

Similar works

Full text

Available Versions

Enlighten