A challenging problem in systems biology is parameter inference in mechanistic models of signalling pathways. In the present article, we investigate an approach based on gradient matching and nonparametric Bayesian modelling with Gaussian processes. We evaluate the method on two biological systems, related to the regulation of PIF4/5 in Arabidopsis thaliana, and the JAK/STAT signal transduction pathway

Cretella, R.

Dondelinger, F.

Filippone, M.

Husmeier, D.

Millar, A.

Palmer, T.

Rogers, S.

Smith, R.

English

Enlighten

Enlighten – Research publications by members of the University of Glasgow http://eprints.gla.ac.uk        Dondelinger, F., Rogers, S. , Filippone, M. , Cretella, R., Palmer, T., Smith, R., Millar, A., and Husmeier, D. (2012) Parameter inference in mechanistic models of cellular regulation and signalling pathways using gradient matching. In: WCSB2012 - 9th International Workshop on Computational Systems Biology, 4-6 Jun 2012, Ulm, Germany.  http://eprints.gla.ac.uk/65591/  Deposited on: 20 June 2012   PARAMETER INFERENCE IN MECHANISTIC MODELS OF CELLULARREGULATION AND SIGNALLING PATHWAYS USING GRADIENT MATCHINGFrank Dondelinger1, Simon Rogers2, Maurizio Filippone2, Roberta Cretella3, Tim M. Palmer3,Robert W. Smith4, Andrew J. Millar4 and Dirk Husmeier11 School of Mathematics and Statistics, University of Glasgow2 School of Computing Science, University of Glasgow3 Institute of Cardiovascular and Medical Sciences, University of GlasgowUniversity of Glasgow, Glasgow G12 8QQ, Scotland4 SynthSys, University of Edinburgh, Edinburgh EH9 3JDCorresponding author: Dirk.Husmeier@glasgow.ac.ukABSTRACTA challenging problem in systems biology is parameterinference in mechanistic models of signalling pathways.In the present article, we investigate an approach based ongradient matching and nonparametric Bayesian modellingwith Gaussian processes. We evaluate the method on twobiological systems, related to the regulation of PIF4/5 inArabidopsis thaliana, and the JAK/STAT signal transduc-tion pathway.1. INTRODUCTIONA central problem in computational systems biology is theformulation of a consistent mechanistic model descriptionof signalling pathways and molecular processes of cel-lular regulation. While there have been several attemptsat describing these processes at a qualitative level, properstatistical inference is a much more challenging problem.The approach is based on minimising the discrepancy be-tween measured data, e.g. related to the abundance pro-files of some molecular components, and their simulatedvalues. This discrepancy is related to some metric, whichcan be shown to be defined by the assumed noise model interms of a maximum likelihood approach to inference. Forinstance, minimising the root mean square deviation be-tween measured and simulated data is equivalent to max-imizing the likelihood under the assumption of white ad-ditive Gaussian noise. For further details, see [1].The practical difficulties with this approach are two-fold. First, the likelihood landscape is typically ruggedand multimodal, which calls for some form of annealingscheme. Second, each parameter adaptation requires anumerical solution of the differential equations (ODEs),which is computationally expensive and hence limits thenumber of maximum likelihood (ML) optimization stepsor Markov chain Monte Carlo (MCMC) sampling stepsthat can be carried out at reasonable computational costs.A potential solution to this problem is the approachof gradient matching. The idea is that rather than aim-ing to explicitly solve the ODEs, we seek to minimise thediscrepancy between the gradients inferred from the slopeof the interpolant and those predicted from the coupledsystem of ODEs. The former is defined by the regres-sion model and depends on some smoothness or regular-ization parameters. The latter is defined by the system ofcoupled ODEs and depends on its parameters, whose de-termination is the ultimate objective of inference. Earlierapproaches pursued a two-step approach, in which firstthe interpolant was inferred, and then the ODE parame-ters were inferred by minimizing the discrepancy betweenthe time derivatives predicted from the ODEs and thosepredicted from the slope of the interpolant [2]. The dis-advantage of this approach is that the result of parameterinference critically hinges on the quality of the interpola-tion scheme, which once completed is kept fixed. A bet-ter approach, first suggested in [3], is to allow for somefeedback mechanism by which the system of ODEs canact back on the interpolation scheme. For instance, thesystem of ODEs might only match the slopes of the inter-polant for implausible or a priori unlikely parameter con-figurations, in which case the interpolant should be ad-justed. Hence, in order to be viable, the mismatch be-tween the time derivative predicted from the slope of theinterpolant and the one obtained from the ODEs should besystematically reduced in an iterative loop, whereby boththe ODEs and the smoothness parameters of the regres-sion are adapted simultaneously. A variation of this ap-proach was presented in [4], where the authors employedparallel tempering using the data-smoothing hyperparam-eter, which improves sampling efficiency.2. METHODOur approach is based on non-parametricBayesian regres-sion with Gaussian processes, following Calderhead et al.in [5]. Their approach draws on the fact that the derivativeof a Gaussian process is a Gaussian process again. Thisrenders Gaussian process regression a natural tool for thedouble objective of nonlinear data interpolation and gradi-ent matching along the lines outlined in the previous sec-tion. The result is a hierarchical Bayesian model whichallows for parameter inference. We will have a systemFigure 1. Outline of the approach of gradient matching. Left panel: the slope of the interpolant informs the inferenceof the ODE parameters θ. Centre panel: the derivative dxdtpredicted based on the updated parameters from the ODEsf(x(t), θ) is fed back to the interpolation scheme. Right panel: by iterating the above two steps, the mismatch betweenthe time derivative predicted from the slope of the interpolant and the one obtained from the ODEs is systematicallyreduced, thereby adapting both the ODE parameters and the smoothness parameters of the regression simultaneously.0 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0 Gap 10 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0.1 Gap 10 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0.2 Gap 10 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0 Gap 10 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0.1 Gap 10 5 10 15 200.00.40.8TimepointsPIF4/5TrueSimulatedSampledNoise 0.2 Gap 1Figure 2. PIF4/5 expression sampled every hour, with varying observational noise. We show the true (noiseless) expres-sion values, versus expression values simulated from the ODE system using mean sampled parameters. We also show thevalues of PIF4/5 sampled from the GP model. Top: Calderhead et al. model [5]. Bottom: Improved model. Note that thedashed line in the top panel is out of scale.of coupled ODEs, which predict the time derivative; theODE parameters are adapted so as to minimize the devia-tion from the time derivatives predicted with the Gaussianprocess. A novel aspect of our approach, which consti-tutes an important improvement on the method proposedin [5], is the fact that the smoothness hyper-parametersof the Gaussian process are adapted simultaneously alongwith the parameters of the ODEs. A mathematical de-scription of this approach is beyond the scope and pagelimit of this paper and will be presented elsewhere. Aschematic representation is depicted in Figure 1.3. DATAWe test our approach on two biological systems; gene reg-ulation in the circadian clock of Arabidopsis thaliana, andreceptor signal transduction in the JAK/STAT pathway.3.1. The PIF4/5 modelWe apply our GP parameter inference method to a modelfor gene regulation of genes PIF4 and PIF5 by TOC1 inthe circadian clock gene regulatory network of Arabidop-sis thaliana. The overall network is represented by theLocke 2-loop model [6], with fixed parameters set fol-lowing [7]. Only the parameters involved in regulationof PIF4 and PIF5 are inferred. As the expression profilesare very similar, we simplify the model to represent genesPIF4 and PIF5 as a combined gene PIF4/5. We are in-terested in the promoter strength s, the rate constant Kdand Hill coefficient h of the regulation by TOC1, and thedegradation rate d of the PIF4/5 mRNA. The regulationis represented by the following ODE:d[PIF4/5]dt= s ·KhdKhd+ [TOC1]− d · [PIF4/5] (1)where [PIF4/5] and [TOC1] represent the concentrationof PIF4/5 and TOC1, respectively, and t represents time.3.2. The JAK/STAT pathwayWe analyse a model for interleukin-6 signalling (IL-6) invascular endothelial cells. IL-6 binds to a receptor on the0 100 200 300 400 500 60001234TimeConcentrationSpecies: RParam. Error: 0.001 0.001 0 0.022Parameter: k1f Species: RParameter ValueDensity0.00 0.05 0.10 0.15 0.20 0.250246812Parameter: k14 Species: RParameter ValueDensity0.0 0.2 0.4 0.6 0.80123450 100 200 300 400 500 6000.00.51.01.52.0TimeConcentrationSpecies: SOCS3.R*Param. Error: 0.007 0.062 0.004Parameter: k11 Species: SOCS3.R*Parameter ValueDensity0.00 0.05 0.10 0.15 0.2002468Parameter: k10b Species: SOCS3.R*Parameter ValueDensity0.00 0.05 0.10 0.15 0.20024680 100 200 300 400 500 60005152535TimeConcentrationSpecies: 2−STAT3*_NParam. Error: 0 0.002 0 0 0.001 0Parameter: k6f Species: 2−STAT3*_NParameter ValueDensity0.000 0.005 0.010 0.015 0.020 0.025 0.030020406080Parameter: k8 Species: 2−STAT3*_NParameter ValueDensity0.00 0.02 0.04 0.06 0.08 0.10 0.12051020300 100 200 300 400 500 6000.01.02.03.0TimeConcentrationSpecies: STAT3.R*Param. Error: 0.005 0.001 0.138 0.001Parameter: k3 Species: STAT3.R*Parameter ValueDensity0.0 0.2 0.4 0.6051015Parameter: k13 Species: STAT3.R*Parameter ValueDensity0.0 0.2 0.4 0.6051015Figure 3. JAK/STAT pathway results. Row 1: R. Row 2: SOCS3.R∗. Row 3: 2STAT 3∗ N . Row 4: STAT 3.R∗. See(2) for the equations. True (solid line) and inferred (points with error bars) concentrations for each species. The inferredvalues were obtained by drawing parameter samples from the posterior and forward-simulating the ODE system withthese parameters. The error bars show one standard deviation. Histograms show the distribution of sampled parameters,the dashed line represents the gamma prior on the parameters in the hierarchical Bayesian model, and the horizontal lineindicates the true value of the parameter.plasma membrane, activating the JAK/STAT pathway [8].The receptor is phosphorylated, creating docking sites forsignalling molecules like STAT3. STAT3 binds to the phos-phorylated receptor and is phosphorylated itself. Phos-phorylated STAT3 molecules are released from the recep-tor, dimerize and then migrate to the nucleus to triggermRNA transcription of target proteins like SOCS3. SOCS3acts as a feedback mechanism for the signalling pathway:it binds to active receptors to prevent STAT3 activation andto provide a signal termination. The model we consideris a complex system comprising 13 species and 19 pa-rameters. The dynamics of the system are described bymass-action kinetics, with non-linear interactions amongspecies. Under the assumption of full observation of allspecies, we can decompose the system into 13 subsys-tems, one per species. This simplifies inference, and al-lows us to investigate the local identifiability of parame-ters in this model. For space reasons, we only reproducethe subset of equations consisting of the species in Fig-ure 3 below; the full system will be presented in a futurepaper.d[R]dt= −kf1[R] + kb1[R∗] + k11[SOCS3.R∗] +k14[SOCS3.STAT3.R∗]d[SOCS3.R∗]dt= −k11[SOCS3.R∗] + kf10[SOCS3][R∗]−kb10[SOCS3.R∗]d[2STAT3∗ N ]dt= −kf6[2STAT3∗ N ][P300] +kb6[2STAT3∗N.P300] +k5[2STAT3∗C]− k7[2STAT3∗N ] +k8[2STAT3∗N.P300]d[STAT3.R∗]dt= kf2[STAT3][R∗]− kb2[STAT3.R∗]−k13[SOCS3][STAT3.R∗]−k3[STAT3.R∗] (2)Here the square brackets, [.], indicate concentrations,and the ki are kinetic parameters.4. RESULTSWe simulated data from the PIF4/5 model using the pa-rameters (s = 1,Kd = 0.46, h = 2, d = 1), addingobservation noise with standard deviation in {0, 0.1, 0.2}.The parameter inference was done by sampling from theposterior using MCMC with the model from [5], as wellas using MCMCwith the improvedmodel from Section 2.The time period of the measurements is 24 hours, and theinterval (gap) between observed points is 1 hour, result-ing in 25 datapoints. Figure 2 shows the noiseless con-centrations sampled from the GP model, and simulatedPIF4/5 concentrations from the true parameters and fromthe sampled parameters, for MCMC with the Calderheadet al. model in [5], and MCMC with the improved model.For the JAK/STAT pathway, we generated simulationdata for 600 timepoints with parameter values that gaverealistic behaviour for the different species. The data wassampled at intervals of 60 seconds from time zero, making11 timepoints in total. Due to space restrictions, we onlypresent results for a subset of species. Figure 3 shows theresults for the inactive receptor R, for the SOCS3/R* com-plex (where R* is the activated receptor), for the activated2-STAT3 complex in the nucleus and for the STAT3/R*complex.5. DISCUSSIONOur results demonstrate that gradientmatching is a promis-ing approach for parameter estimation in ODE systems.We have demonstrated that our improvement on themethoddescribed in [5] produces better results in the PIF4/5 sys-tem (Figure 2) in the presence of observation noise.Our application to the JAK/STAT signalling pathwayshows that the approach we have taken is promising, butthat some challenges remain. It also allows us to drawsome inferences about the properties of this system. Forsome species, such as R and SOCS3.R* in rows 1 and 2 ofFigure 3, we obtain very good predictions for the concen-trations of the species, as well as giving a good estimatefor some of the parameters, such as k1f and k11. However,we can see that there is a lot of uncertainty about the in-ferred parameters, even though the species concentrationsare predicted quite well. This points to a problem withlack of identifiability in parameter space, related to ridgesin the likelihood. For example, if the rate limiting chem-ical kinetics depend on the ratio of two kinetic constants,then the confidence or credible intervals of the individualparameters may be large without that being reflected bythe prediction uncertainty, as long as the posterior distri-bution of the ratio is peaked.We notice that for species 2-STAT3* N (row 3) in Fig-ure 3, the posterior probability mass for the parameters isin the tail of the prior distribution. This implies two things.First, the data are informative with respect to the infer-ence of some parameters. Second, the prior has not beenchosen very well for this example, and should be chosenmore appropriately. Finally, the bottom row of Figure 3shows a mismatch between the true and predicted signalof STAT3.R*. This points to an intrinsic difficulty withthe gradient matching approach. The mismatch is mainlydue to a short transient region around time point 100. Inthe following time segment the gradient is well matched,whereas the signal itself shows a strong deviation. This in-dicates that in scenarios of this form, the likelihood land-scape for the explicit solution of the ODEs differs system-atically from the one obtained with gradient matching.In conclusion, our gradientmatching approach demon-strates good predictions for realistic biological systems,and promises to become a useful tool for parameter esti-mation in system biology. However, some challenges re-lated to unidentifiable parameters and non-stationary sig-nals still remain.6. ACKNOWLEDGEMENTSThis work was funded by Bridging-the-Gap EPSRC grant59229/1. R.W.S. was supported by BBSRC (BB/F59011/1,BB/F005237/1). SynthSys is a Centre for Integrative andSystems Biology partly supported by BBSRC and EPSRC(BB/G019621).7. REFERENCES[1] N. D. Lawrence, M. Girolami, M. Rattray, andG. Sanguinetti, Eds., Learning and Inference in Com-putational Systems Biology, MIT Press, 2010.[2] A.A. Poyton, M.S. Varziri, K.B. McAuley, P.J.McLellan, and J.O. Ramsay, “Parameter estimationin continuous-time dynamic models using principaldifferential analysis,” Computers & chemical engi-neering, vol. 30, no. 4, pp. 698–708, 2006.[3] J.O. Ramsay, G. Hooker, D. Campbell, and J. Cao,“Parameter estimation for differential equations: ageneralized smoothing approach,” Journal ofthe Royal Statistical Society: Series B (StatisticalMethodology), vol. 69, no. 5, pp. 741–796, 2007.[4] D. Campbell and R.J. Steele, “Smooth functionaltempering for nonlinear differential equation models,”Statistics and Computing, pp. 1–15, 2012.[5] B. Calderhead, M. Girolami, and N.D. Lawrence,“Accelerating Bayesian inference over nonlinear dif-ferential equations with Gaussian processes,” Neu-ral Information Processing Systems (NIPS), vol. 22,2008.[6] J.C.W. Locke, M.M. Southern, L. Kozma-Bognar,V. Hibberd, P.E. Brown, M.S. Turner, and A.J. Millar,“Extension of a genetic network model by iterativeexperimentation and mathematical analysis,” Molec-ular Systems Biology, vol. 1, pp. (online), 2005.[7] A. Pokhilko, S.K. Hodge, K. Stratford, K. Knox, K.D.Edwards, A.W. Thomson, T. Mizuno, and A.J. Mil-lar, “Data assimilation constrains new connectionsand components in a complex, eukaryotic circadianclock model,” Molecular systems biology, vol. 6, no.1, 2010.[8] P.C. Heinrich, I. Behrmann, G. Mu¨ller-Newen,F. Schaper, and L. Graeve, “Interleukin-6-type cy-tokine signalling through the gp130/Jak/STAT path-way.,” Biochemical Journal, vol. 334, no. Pt 2, pp.297, 1998.

Accelerating Bayesian inference over nonlinear differential equations with Gaussian processes,”

Data assimilation constrains new connections and components in a complex, eukaryotic circadian clock model,” Molecular systems biology,

Extension of a genetic network model by iterative experimentation and mathematical analysis,”

Interleukin-6-type cytokine signalling through the gp130/Jak/STAT pathway.,”

Parameter estimation for differential equations: a generalized smoothing approach,”

Parameter estimation in continuous-time dynamic models using principal differential analysis,”

Parameter inference in mechanistic models of cellular regulation and signalling pathways using gradient matching

Enlighten: Publications

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

Parameter inference in mechanistic models of cellular regulation and signalling pathways using gradient matching

Abstract

Similar works

Full text

Available Versions

Enlighten

Enlighten: Publications