 # Fast Gaussian Process Based Gradient Matching for Parameter Identification in Systems of Nonlinear ODEs

Parameter identification and comparison of dynamical systems is a challenging task in many fields. Bayesian approaches based on Gaussian process regression over time-series data have been successfully applied to infer the parameters of a dynamical system without explicitly solving it. While the benefits in computational cost are well established, a rigorous mathematical framework has been missing. We offer a novel interpretation which leads to a better understanding and improvements in state-of-the-art performance in terms of accuracy for nonlinear dynamical systems.

## Code Repositories

Sample code for the NIPS paper "Scalable Variational Inference for Dynamical Systems"

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The underlying mechanism of many processes in science and engineering can often be described by ordinary differential equations (ODE). While the form of dynamical systems, the ODEs, can often be derived using expert knowledge, the parameters are usually unknown, can not be directly measured and have to be estimated from empirical time series data. Since nonlinear ODEs typically do not have a closed form solution, standard methods for statistical inference require the computationally expensive numerical integration of the ODEs every time the parameters are changed

To circumvent the high computational cost of numerical integration, gradient matching techniques have been proposed (e.g. Ramsay et al., 2007; Dondelinger et al., 2013; Niu et al., 2016; Gorbach et al., 2017a)

. Gradient matching is based on minimizing the difference between a model interpolating the dynamics of the state variables and the time derivatives provided by the ODEs. The first steps of this approach go back to spline-based methods, with an overview given by

Ramsay et al. (2007). Calderhead et al. (2008) then proposed a fully probabilistic treatment by using Gaussian process models, increasing data efficiency and providing final uncertainty estimates. To match the gradients of the GP model to the gradients provided by the ODEs, Calderhead et al. (2008)

introduce a product of experts heuristics (PoE). This heuristic has since become the state-of-the-art explanation, being reused, e.g., in

Dondelinger et al. (2013) and Gorbach et al. (2017a). The advantages of the method have been well established, with applications in biomedical domains e.g., systems biology or neuroscience (Babtie et al., 2014; Macdonald and Husmeier, 2015; Pfister et al., 2018). Given these applications, a solid understanding and theoretical framework of theses methods is of critical importance.

However, there are two important open issues. Regarding the theoretical framework, there has been some criticism regarding the use of the PoE in this context, e.g., by Wang and Barber (2014), who introduced a different modeling paradigm. Regarding empirical performance, Gorbach et al. (2017a) recently introduced a mean field variational scheme to decrease the run time. However, despite the expected run time savings, Gorbach et al. (2017a) also reported a significantly increased accuracy of the variational approach compared to the MCMC sampling based counterparts.

While the criticism of the product of experts approach (Wang and Barber, 2014) lead to the controversy of mechanistic modeling with Gaussian processes (Macdonald et al., 2015) where the theoretical inconsistencies of the modeling proposal of Wang and Barber (2014) have been outlined, a similar analysis investigating the theoretical foundation of the state of the art approach underlying all current modeling proposals (Calderhead et al., 2008; Dondelinger et al., 2013; Gorbach et al., 2017a) has been missing.

Our contributions. In this work, we

1. analyze the Product of Experts (PoE) heuristic, discovering and explaining theoretical inconsistencies of the state of the art approach,

2. provide a theoretical framework, replacing the criticized product of experts heuristic,

3. identify the cause of the performance gains of the variational approach of Gorbach et al. (2017a) over sampling-based methods,

4. combine these insights to create a novel algorithm improving on state-of-the-art performance for nonlinear systems in terms of accuracy and robustness, while providing a more computationally efficient sampling scheme reducing its run time by roughly 35%.111Code publicly available at

## 2 Preliminaries

### 2.1 Problem Formulation

In this work, we analyze an arbitrary system of parametric ordinary differential equations, whose state derivatives can be parametrized by a time independent parameter vector

. In this setting, the true evolution of the dynamical system is given by

 ˙x=f(x,θ) (1)

where are the time dependent states of the system and is an arbitrary, potentially nonlinear vector valued function.

While Equation (1) is meant to represent the dynamics of the system at all time points, , and will be used throughout this paper to denote the vector containing the time evolution of one state or one state derivative at the observation times , i.e. etc.

At N discrete time points , the state trajectories are observed with additive, i.i.d. Gaussian noise , i.e.,

 y(ti)=x(ti)+ϵ(ti),i=1…N, (2)

or equivalently

 p(y|x,σ)=N(y|x,σ2I). (3)

Given these noisy observations and the functional form of , the goal is to infer the true parameters .

For the sake of clarity, we present the theory for a one-dimensional system only and assume that all observations were created using one realization of the experiment. The extension to multidimensional systems (as done in the experiments, Section 5) or repeated experiments is straightforward but omitted to simplify the notation.

### 2.2 Modeling

The key idea of GP-based gradient matching is to build a GP regression model mapping the time points to the corresponding state values. For this, one needs to choose an appropriate kernel function

, which is parametrized by the hyperparameters

. Both, the choice of kernels as well as how to fit its hyperparameters is discussed in Section 7.3.

Once the kernel and its hyperparameters are fixed, the covariance matrix , whose elements are given by , can be constructed and used to define a standard zero mean Gaussian process prior on the states:

 p(x|ϕ)=N(x|0,Cϕ). (4)

As differentiation is a linear operation, the derivative of a Gaussian process is again a Gaussian process. Using probabilistic calculus (see supplementary material, Section 7.1 for details), this fact leads to a distribution over the derivatives conditioned on the states at the observation points:

 p(˙x|x,ϕ)=N(˙x|Dx,A). (5)

Lastly, the information provided by the differential equation is used as well. For known states and parameters, one can calculate the derivatives using equation (1

). A potential modeling mismatch between the output of the ODEs and the derivatives of the GP model is accounted for by introducing isotropic Gaussian noise with standard deviation

, leading to the following Gaussian distribution over the derivatives:

 p(˙x|x,θ,γ)=N(˙x|f(x,θ),γI). (6)

The modeling assumptions are summarized in the graphical models shown in Figure 1.

### 2.3 Inference

As stated in Section 2.1, the main goal of the inference process is to learn the parameters using the noisy observations . Thus, it is necessary to connect the two graphical models shown in Figure 1. As shown in the previous section, and the

represent the same variables in both models. However, it is not straightforward to use this fact to combine the two. While it is intuitive to use the probability density over

of the Gaussian process model directly as the prior for in the ODE response model, handling is more challenging. In both models, is a dependent variable. Thus, some heuristic is needed to combine the two conditional distributions and .

#### 2.3.1 Product of experts heuristic

The main idea of the product of experts, originally introduced by Hinton (2002), is to infer the probability density of a variable by normalizing the product of multiple expert densities. Calderhead et al. (2008) use this to connect the two distributions over , leading to

 p(˙x|x,ϕ,θ,γ)∝p(˙x|x,ϕ)p(˙x|x,θ,γ) (7)

The idea of this approach is that the resulting density only assigns high probability if both experts assign high probabilities. Hence, it considers only cases in which both experts agree. It is thus based on the intuition that the true should correspond to a model that agrees both with the ODE model and the observed data. While this is intuitively well-motivated, we will show that the product of experts heuristic leads to theoretical difficulties and offer an alternative.

#### 2.3.2 Markov Chain Monte Carlo based methods

Calderhead et al. (2008) combine the product of experts with equations (3), (4) and (5) and some suitable prior over

to obtain a joint distribution

. After integrating out , which can be done analytically since Gaussian processes are closed under linear operators (and using some proportionality arguments), a sampling scheme was derived that consists of two MCMC steps. First, the hyperparameters of the GP, and , are sampled from the conditional distribution . Then, a second MCMC scheme is deployed to infer the parameters of the ODE model, and , by sampling from the conditional distribution .

Dondelinger et al. (2013) then reformulated the approach by directly calculating the joint distribution

 p(y,x,θ,ϕ,γ,σ) ∝p(θ) ×N(x|0,Cϕ) ×N(y|x,σ2I) ×N(f(x,θ)|Dx,A+γI), (8)

where the proportionality is meant to be taken w.r.t. the latent states and the ODE parameters . Here denotes some prior on the ODE parameters. This approach was named Adaptive Gradient Matching ().

#### 2.3.3 Variational inference

The main idea of Variational Gradient Matching (), introduced by Gorbach et al. (2017a), is to substitute the MCMC inference scheme of with a mean field variational inference approach, approximating the density in Equation (8) with a fully factorized Gaussian over the states and the parameters . To obtain analytical solutions, the functional form of the ODEs is restricted to locally linear functions that could be written as

 f(x,θ)=∑iθi∏j∈MixjwhereMi⊆{1,...,K}. (9)

As perhaps expected, is magnitudes faster than the previous sampling approaches. However, despite being a variational approach, was also able to provide significantly more accurate parameter estimates than both sampling-based approaches of Calderhead et al. (2008) and Dondelinger et al. (2013). In Section 4, we provide justification for these surprising performance differences.

## 3 Theory

### 3.1 Analysis Of The Product Of Experts Approach

As previously stated, Calderhead et al. (2008), Dondelinger et al. (2013) and Gorbach et al. (2017a) all use a product of experts to obtain as stated in Equation (7).

In this section, we will first provide an argument based on graphical models and then an argument based on the original mathematical derivation to illustrate challenges arising from this heuristic.

Figure 2 depicts what is happening if the product of experts approach is applied in the gradient matching framework. Figure 1(a) depicts the graphical model after the two models have been merged using the product of experts heuristic of Equation (7). Using the distribution over of the Gaussian process model 0(a) as a prior for the in the ODE response model 0(b), effectively leads to merging the two nodes representing . Furthermore, the product of experts heuristic implies by its definition that after applying Equation (7), is only depending on , , and .

In the graphical model in Figure 1(a), the problem is already visible. The ultimate goal of merging the two graphical models is to create a probabilistic link between the observations and the ODE parameters . However, the newly created connection between these two variables is given via , which has no outgoing edges and of which no observations are available. Marginalizing out as proposed in the traditional approaches consequently leads to the graphical model in Figure 1(b). As there is no directed path connecting other variables via , all the different components are now independent. Consequently, the posterior over is now given by the prior we put on in the first place.

This problem can further be illustrated by the mathematical derivations in the original paper of Calderhead et al. (2008). After calculating the necessary normalization constants, the last equation in the third chapter is equivalent to stating

 p(θ,γ|x,ϕ,σ)=∫p(θ)p(γ)p(˙x|x,θ,γ,ϕ,σ)d˙x. (10)

It is clear that this equation should simplify to

 p(θ,γ|x,ϕ,σ)=p(θ)p(γ). (11)

Thus, one could argue that any links that are not present in the graphical model of Figure 1(b) but found by Calderhead et al. (2008) and reused in Dondelinger et al. (2013) and Gorbach et al. (2017a) were created by improper normalization of the density .

### 3.2 Adapting The Original Graphical Model

Despite these technical difficulties arising from the PoE heuristic, the approaches provide good empirical results and have been used in practice, e.g., by Babtie et al. (2014). In what follows, we derive an alternative model and mathematical justification for Equation (8) to provide a theoretical framework explaining the good empirical performance of Gaussian process-based gradient matching approaches, especially from Gorbach et al. (2017a), which uses only weak or nonexistent priors.

The graphical model shown in Figure 3 offers an alternative approach to the product of experts heuristic. The top two layers are equivalent to a GP prior on the states, the induced GP on the derivatives and the observation model, as shown in Figure 0(a).

The interesting part of the new graphical model is the bottom layer. Instead of adding a second graphical model like in Figure 0(b)

to account for the ODE response, two additional random variables are introduced.

is the deterministic output of the ODEs, assuming the values of and are given, i.e. . The deterministic nature of this equation is represented as a Dirac delta function:

 (12)

If the GP model were able to capture the true states and true derivatives perfectly, this new random variable should be equivalent to the derivatives of the GP, i.e., . However, to compensate for a potential model mismatch and slight errors of both GP states and GP derivatives, this condition is relaxed to

 F1=˙x+ϵ=:F2, ϵ∼N(0,γI) (13)

In the graphical model, this intuitive argument is encoded via the random variable . Given the provided by the GP model, Gaussian noise with standard deviation is added to create , whose probability density can thus be described as

 p(F2|˙x,γ)=N(F2|˙x,γI). (14)

The equality constraint given by equation (13) is represented in the graphical model by the undirected edge between and . When doing inference, this undirected edge is incorporated in the joint density via a Dirac delta function . Thus, the joint density of the graphical model represented in Figure 3 can be written as

 p(x,˙x,y, F1,F2,θ|ϕ,σ,γ)=p(θ) ×p(x|ϕ)p(˙x|x,ϕ)p(y|x,σ) ×p(F1|θ,x)p(F2|˙x,γI)δ(F1−F2). (15)

### 3.3 Inference In The New Model

Given all the definitions in the previous section, inference can now be directly performed without the need for additional heuristics. The result is a theoretically sound justification of the main result of Calderhead et al. (2008) and Dondelinger et al. (2013):

###### Theorem

Given the modeling assumptions summarized in the graphical model in Figure 3,

 p(x,θ|y,ϕ,γ,σ) ∝p(θ) ×N(x|0,Cϕ) ×N(y|x,σ2I) ×N(f(x,θ)|Dx,A+γI). (16)

The proof can be found in the supplementary material, section 7.2.

## 4 Fast Gaussian Process Gradient Matching

### 4.1 Hyperparameters

Using the theoretical framework of the previous section, it is clear that the performance of any algorithm based on equation 16 will heavily rely on the quality of the hyperparameters , as , and are all depending on . In GP regression, it is common to fit the hyperparameters to the observations using a maximum likelihood scheme (see supplementary material, section (7.3)). However, neither nor actually do this.

In , the hyperparameters are inferred concurrently to the states and the ODE parameters

in one big MCMC scheme. Besides the need for clear hyperpriors on the hyperparameters and obvious drawbacks regarding running time, e.g. at each MCMC step

, and have to be recalculated and potentially inverted, this leads to a rough probability landscape requiring a complicated multi-chain setup (Dondelinger et al., 2013).

In , this problem was sidestepped by treating the hyperparameters as tuning parameters that had to be set by hand. As these hyperparameters are not learned from data, it is obviously not a fair comparison. In our experiments, we thus used a maximum likelihood scheme to infer the hyperparameters before using the implementation of Gorbach et al. (2017a). To indicate this modification, the modified VGM algorithm was called .

### 4.2 Fgpgm

However, is still outperforming significantly as shown in the experiments in section 5. This suggests that the concurrent optimization of , and suggested by Dondelinger et al. (2013) is not helpful for the performance of the system. Based on this insight, we propose the sequential approach shown in Algorithm 1. In a first step, the Gaussian process model is fit to the standardized data by calculating the hyperparameters via equation (7.3.1). Then, the states and ODE parameters are inferred using a one chain MCMC scheme on the density given by equation (8).

## 5 Experiments

For all experiments involving , the R toolbox deGradInfer (Macdonald and Dondelinger, 2017) published alongside (Macdonald, 2017) was used. The toolbox was provided with code to run two experiments, namely Lotka Volterra and Protein Transduction. Both of these systems are used in this paper, as they are the two standard benchmark systems used in all previous publications. It should be noted however that the applicability of is not restricted to these systems. Unlike , refrains from using hard to motivate hyperpriors and our publicly available implementation can easily be adapted to new settings.

For details about implementation and additional plots, refer to the supplementary material in section 7.6. It should be noted however that all algorithms were provided with one hyperparameter. While the toolbox of had to be provided with the true standard deviation of the observation noise, and were provided with . The was determined by testing eight different values logarithmically spaced between 1 and and comparing the results based on observation fit. For more details regarding the experimental setup, we refer to the appendix of Wenk et al. (2019).

### 5.1 Lotka Volterra

The first system in question is the Lotka Volterra system originally proposed in Lotka (1978). It describes a two dimensional system whose dynamics are given by

 ˙x1(t) =θ1x1(t)−θ2x1(t)x2(t) ˙x2(t) =−θ3x2(t)+θ4x1(t)x2(t)

We reproduce the experimental setup of Gorbach et al. (2017a), i.e., the system was observed in the time interval at 20 evenly spaced observation times. The true states were initialized with and Gaussian noise with standard deviation 0.1 (low noise) and standard deviation 0.5 (high noise) was added to create the observations.

### 5.2 Protein Transduction

The second system to be analyzed is called Protein Transduction and was originally proposed in Vyshemirsky and Girolami (2008). It is known to be notoriously difficult to fit with unidentifiable parameters (Dondelinger et al., 2013). The dynamics are given by:

 ˙S =−θ1S−θ2SR+θ3RS ˙dS =θ1S ˙R =−θ2SR+θ3RS+θ5Rppθ6+Rpp ˙RS =θ2SR−θ3RS−θ4RS ˙Rpp =θ4RS−θ5Rppθ6+Rpp (17)

It should be noted that these dynamics contain nonlinear terms violating the functional form assumption given in equation (9) of inherited by . Nevertheless, both and can still be applied. For , was set to , while was provided with the true observation noise standard deviations. The experimental setup of Dondelinger et al. (2013) was copied, i.e. the system was observed in the time interval at the discrete observation times , the states were initialized with and the parameters were set to . As in Dondelinger et al. (2013), Gaussian noise was added with standard deviation 0.001 (low noise) and 0.01 (high noise). As in the previous papers, a sigmoid kernel was used to deal with the logarithmically spaced observation times and the typically spiky form of the dynamics. Figure 5: States after numerical integration of the inferred parameters in the high noise case of Lotka Volterra. Ground truth (red), median (black) and 75% quantiles (gray) over 100 independent noise realizations.

### 5.3 Evaluation

As the parameters of the Protein Transduction system are unidentifiable, comparing parameter values is not a good metric to rank approaches. Instead, after running the algorithms, we used a numerical integrator to obtain the trajectories corresponding to the dynamical system whose dynamics is given by the the inferred parameters. Then, the RMSE of these trajectories compared to the ground truth at the observation times were evaluated. For each experimental setting, this procedure was repeated for 100 different noise realizations. The results are shown in Figure 4. In Figure 5, we show median plots for the high noise setting of Lotka Volterra, while the running time between the state of the art and is shown in Figure 6.

While achieving run time savings of 35%, shows an increased accuracy, reducing the state RMSE up to 62 % as shown in Table 1

and shows much lower variability across noise realizations. The effect is especially striking for the nonlinear system Protein Transduction, where the variance of the states

and has been reduced by at least one order of magnitude, while reducing the bias by more than 30%, see Figures 3(c) and 3(d). (a) run time LV Figure 7: One input sample and median plots of the numerically integrated states after parameter inference for the FHN system with SNR 10. Ground truth (red), median (black) and 75% quantiles (gray) over 100 independent noise realizations. In the most left plots, the black dots represent the noisy observations. Figure 8: Median plots of the first state of FHN with SNR 100 in case of 10, 25, 50 and 100 observations. Same colors as in Figure 7. Plots of state 2 can be found in the supplementary material.

### 5.4 Spiky Dynamics

As can be seen in Figure 5, all GP based gradient matching algorithms converge to parameter settings where the trajectories are smoother than the ground truth. While learning the hyperparameters in a pre-processing step clearly reduces this effect for and , there is still some bias. If only few observations are available, the GP prior on the states tends to smooth out part of the system dynamics. This ”smoothing bias” is then passed on to the ODE parameters in the gradient matching scheme.

To investigate the practical importance of this effect, we evaluate on a third system proposed by FitzHugh (1961) and Nagumo et al. (1962)

for modeling giant squid neurons. Abbreviating the name of its inventors, we will refer to it as the FHN system. Its dynamics are given by the two ODEs

 ˙V =θ1(V−V33+R) (18) ˙R =1θ1(V−θ2+θ3R) (19)

Due to its highly nonlinear terms, the FHN system has notoriously fast changing dynamics. One example realization including noisy observations is shown in Figure 7. To account for the spikier behavior of the system, we used a Matérn52 kernel. was set to .

While the bias towards smoother trajectories is clearly visible, finds parameters that tightly hug the ground truth. As to be expected, the smoothing bias gets smaller if we add more observations. Furthermore, increasing the SNR to 100 as shown in Figure 8 leads to an impressive accuracy, even if we reduce the amount of observations to just 10, especially as is a hyper-prior free, statistical method

## 6 Discussion

Gradient matching is a successful tool to circumvent the computational cost of numerical integration for Bayesian parameter identification in dynamical systems, especially if the dynamics, like in most real world systems, are reasonably smooth. Previous Gaussian process-based approaches used a criticized product of experts heuristics, which leads to technical difficulties. We illustrated these theoretical problems and provided a novel, sound formulation that does not rely on a PoE. We furthermore explained the surprising performance gains of variational over sampling-based approaches and then combined these insights to propose a new algorithm called , which jointly learns states and parameters with improved state-of-the-art performance in terms of accuracy, run time, robustness and reduction of ”smoothing bias” for general nonlinear dynamical systems.

Unlike the previous MCMC approaches, uses a one-chain Metropolis Hastings scheme, which is much easier to tune than the previously used, complicated multi chain setups. Due to the sequential fitting of the hyperparameters, the wild behavior of the probability density motivating the setup in Dondelinger et al. (2013, section 3) and the need for hard to motivate hyperpriors is successfully avoided. Consequently, the inference is significantly more accurate, faster, more robust and significantly decreases the smoothing bias of GP based gradient matching schemes.

#### Acknowledgements

This research was partially supported by the Max Planck ETH Center for Learning Systems and SNSF grant 200020_159557.

## References

• Babtie et al. (2014) Ann C Babtie, Paul Kirk, and Michael PH Stumpf. Topological sensitivity analysis for systems biology. Proceedings of the National Academy of Sciences, 111(52):18507–18512, 2014.
• Bauer et al. (2017) Stefan Bauer, Nico S Gorbach, Djordje Miladinovic, and Joachim M. Buhmann. Efficient and flexible inference for stochastic systems. Neural Information Processing Systems (NIPS), 2017.
• Bishop (2006) Christopher M Bishop. Springer-Verlag New York, 2006.
• Calderhead et al. (2008) Ben Calderhead, Mark Girolami  and Neil D. Lawrence.

Accelerating bayesian inference over nonliner differential equations with gaussian processes.

Neural Information Processing Systems (NIPS), 2008.
• Choi and Schervish (2007) Taeryon Choi and Mark J Schervish. On posterior consistency in nonparametric regression problems.

Journal of Multivariate Analysis

, 98(10):1969–1987, 2007.
• Dondelinger et al. (2013) Frank Dondelinger, Maurizio Filippone, Simon Rogers  and Dirk Husmeier. Ode parameter inference using adaptive gradient matching with gaussian processes.

International Conference on Artificial Intelligence and Statistics (AISTATS)

, 2013.
• Duvenaud (2014) David Duvenaud. Automatic model construction with Gaussian processes. PhD thesis, University of Cambridge, 2014.
• Duvenaud et al. (2013) David Duvenaud, James Robert Lloyd, Roger Grosse, Joshua B Tenenbaum, and Zoubin Ghahramani. Structure discovery in nonparametric regression through compositional kernel search. arXiv preprint arXiv:1302.4922, 2013.
• FitzHugh (1961) Richard FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical journal, 1(6):445–466, 1961.
• González et al. (2014) Javier González, Ivan Vujačić, and Ernst Wit. Reproducing kernel hilbert space based estimation of systems of ordinary differential equations. Pattern Recognition Letters, 45:26–32, 2014.
• Gorbach et al. (2017a) Nico S Gorbach, Stefan Bauer, and Joachim M. Buhmann. Scalable variational inference for dynamical systems. Neural Information Processing Systems (NIPS), 2017a.
• Gorbach et al. (2017b) Nico S Gorbach, Andrew An Bian, Benjamin Fischer, Stefan Bauer, and Joachim M Buhmann. Model selection for gaussian process regression. In German Conference on Pattern Recognition, pages 306–318. Springer, 2017b.
• Hinton (2002) Geoffrey E Hinton.

Training products of experts by minimizing contrastive divergence.

Neural computation, 14(8):1771–1800, 2002.
• Lotka (1978) Alfred J Lotka. The growth of mixed populations: two species competing for a common food supply. In The Golden Age of Theoretical Ecology: 1923–1940, pages 274–286. Springer, 1978.
• Macdonald (2017) Benn Macdonald. Statistical inference for ordinary differential equations using gradient matching. PhD thesis, University of Glasgow, 2017.
• Macdonald and Dondelinger (2017) Benn Macdonald and Frank Dondelinger. deGradInfer: Parameter inference for systems of differential equation.
• Macdonald and Husmeier (2015) Benn Macdonald and Dirk Husmeier. Gradient matching methods for computational inference in mechanistic models for systems biology: a review and comparative analysis. Frontiers in bioengineering and biotechnology, 3, 2015.
• Macdonald et al. (2015) Benn Macdonald, Catherine F. Higham  and Dirk Husmeier. Controversy in mechanistic modemodel with gaussian processes. International Conference on Machine Learning (ICML), 2015.
• Nagumo et al. (1962) Jinichi Nagumo, Suguru Arimoto, and Shuji Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50(10):2061–2070, 1962.
• Niu et al. (2016) Mu Niu, Simon Rogers, Maurizio Filippone, and Dirk Husmeier. Fast inference in nonlinear dynamical systems using gradient matching. International Conference on Machine Learning (ICML), 2016.
• Papoulis and Pillai (2002) Athanasios Papoulis and S Unnikrishna Pillai. Probability, random variables, and stochastic processes. Tata McGraw-Hill Education, 2002.
• Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. Journal of machine learning research, 12(Oct):2825–2830, 2011.
• Petersen et al. (2008) Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008.
• Pfister et al. (2018) Niklas Pfister, Stefan Bauer, and Jonas Peters. Identifying causal structure in large-scale kinetic systems. arXiv preprint arXiv:1810.11776, 2018.
• Ramsay et al. (2007) Jim O Ramsay, Giles Hooker, David Campbell, and Jiguo Cao. Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):741–796, 2007.
• Rasmussen (2004) Carl Edward Rasmussen. Gaussian processes in machine learning. In Advanced lectures on machine learning, pages 63–71. Springer, 2004.
• Vyshemirsky and Girolami (2008) Vladislav Vyshemirsky and Mark A Girolami. Bayesian ranking of biochemical system models. Bioinformatics, 24(6):833–839, 2008.
• Wang and Barber (2014) Yali Wang and David Barber. Gaussian processes for bayesian estimation in ordinary differential equations. International Conference on Machine Learning (ICML), 2014.
• Wenk et al. (2019) Philippe Wenk, Gabriele Abbati, Stefan Bauer, Michael A Osborne, Andreas Krause, and Bernhard Schölkopf. Odin: Ode-informed regression for parameter and state inference in time-continuous dynamical systems. arXiv preprint arXiv:1902.06278, 2019.
• Zhang et al. (2016) Xinxin Zhang, Ruirui Ji, and Yingmin Yi. Inference of transcriptional regulation using multi-objective optimization with gaussian process. In Image and Signal Processing, BioMedical Engineering and Informatics (CISP-BMEI), International Congress on, pages 1820–1825. IEEE, 2016.

## References

• Babtie et al. (2014) Ann C Babtie, Paul Kirk, and Michael PH Stumpf. Topological sensitivity analysis for systems biology. Proceedings of the National Academy of Sciences, 111(52):18507–18512, 2014.
• Bauer et al. (2017) Stefan Bauer, Nico S Gorbach, Djordje Miladinovic, and Joachim M. Buhmann. Efficient and flexible inference for stochastic systems. Neural Information Processing Systems (NIPS), 2017.
• Bishop (2006) Christopher M Bishop. Springer-Verlag New York, 2006.
• Calderhead et al. (2008) Ben Calderhead, Mark Girolami  and Neil D. Lawrence.

Accelerating bayesian inference over nonliner differential equations with gaussian processes.

Neural Information Processing Systems (NIPS), 2008.
• Choi and Schervish (2007) Taeryon Choi and Mark J Schervish. On posterior consistency in nonparametric regression problems.

Journal of Multivariate Analysis

, 98(10):1969–1987, 2007.
• Dondelinger et al. (2013) Frank Dondelinger, Maurizio Filippone, Simon Rogers  and Dirk Husmeier. Ode parameter inference using adaptive gradient matching with gaussian processes.

International Conference on Artificial Intelligence and Statistics (AISTATS)

, 2013.
• Duvenaud (2014) David Duvenaud. Automatic model construction with Gaussian processes. PhD thesis, University of Cambridge, 2014.
• Duvenaud et al. (2013) David Duvenaud, James Robert Lloyd, Roger Grosse, Joshua B Tenenbaum, and Zoubin Ghahramani. Structure discovery in nonparametric regression through compositional kernel search. arXiv preprint arXiv:1302.4922, 2013.
• FitzHugh (1961) Richard FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical journal, 1(6):445–466, 1961.
• González et al. (2014) Javier González, Ivan Vujačić, and Ernst Wit. Reproducing kernel hilbert space based estimation of systems of ordinary differential equations. Pattern Recognition Letters, 45:26–32, 2014.
• Gorbach et al. (2017a) Nico S Gorbach, Stefan Bauer, and Joachim M. Buhmann. Scalable variational inference for dynamical systems. Neural Information Processing Systems (NIPS), 2017a.
• Gorbach et al. (2017b) Nico S Gorbach, Andrew An Bian, Benjamin Fischer, Stefan Bauer, and Joachim M Buhmann. Model selection for gaussian process regression. In German Conference on Pattern Recognition, pages 306–318. Springer, 2017b.
• Hinton (2002) Geoffrey E Hinton.

Training products of experts by minimizing contrastive divergence.

Neural computation, 14(8):1771–1800, 2002.
• Lotka (1978) Alfred J Lotka. The growth of mixed populations: two species competing for a common food supply. In The Golden Age of Theoretical Ecology: 1923–1940, pages 274–286. Springer, 1978.
• Macdonald (2017) Benn Macdonald. Statistical inference for ordinary differential equations using gradient matching. PhD thesis, University of Glasgow, 2017.
• Macdonald and Dondelinger (2017) Benn Macdonald and Frank Dondelinger. deGradInfer: Parameter inference for systems of differential equation.
• Macdonald and Husmeier (2015) Benn Macdonald and Dirk Husmeier. Gradient matching methods for computational inference in mechanistic models for systems biology: a review and comparative analysis. Frontiers in bioengineering and biotechnology, 3, 2015.
• Macdonald et al. (2015) Benn Macdonald, Catherine F. Higham  and Dirk Husmeier. Controversy in mechanistic modemodel with gaussian processes. International Conference on Machine Learning (ICML), 2015.
• Nagumo et al. (1962) Jinichi Nagumo, Suguru Arimoto, and Shuji Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50(10):2061–2070, 1962.
• Niu et al. (2016) Mu Niu, Simon Rogers, Maurizio Filippone, and Dirk Husmeier. Fast inference in nonlinear dynamical systems using gradient matching. International Conference on Machine Learning (ICML), 2016.
• Papoulis and Pillai (2002) Athanasios Papoulis and S Unnikrishna Pillai. Probability, random variables, and stochastic processes. Tata McGraw-Hill Education, 2002.
• Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. Journal of machine learning research, 12(Oct):2825–2830, 2011.
• Petersen et al. (2008) Kaare Brandt Petersen, Michael Syskind Pedersen, et al. The matrix cookbook. Technical University of Denmark, 7(15):510, 2008.
• Pfister et al. (2018) Niklas Pfister, Stefan Bauer, and Jonas Peters. Identifying causal structure in large-scale kinetic systems. arXiv preprint arXiv:1810.11776, 2018.
• Ramsay et al. (2007) Jim O Ramsay, Giles Hooker, David Campbell, and Jiguo Cao. Parameter estimation for differential equations: a generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):741–796, 2007.
• Rasmussen (2004) Carl Edward Rasmussen. Gaussian processes in machine learning. In Advanced lectures on machine learning, pages 63–71. Springer, 2004.
• Vyshemirsky and Girolami (2008) Vladislav Vyshemirsky and Mark A Girolami. Bayesian ranking of biochemical system models. Bioinformatics, 24(6):833–839, 2008.
• Wang and Barber (2014) Yali Wang and David Barber. Gaussian processes for bayesian estimation in ordinary differential equations. International Conference on Machine Learning (ICML), 2014.
• Wenk et al. (2019) Philippe Wenk, Gabriele Abbati, Stefan Bauer, Michael A Osborne, Andreas Krause, and Bernhard Schölkopf. Odin: Ode-informed regression for parameter and state inference in time-continuous dynamical systems. arXiv preprint arXiv:1902.06278, 2019.
• Zhang et al. (2016) Xinxin Zhang, Ruirui Ji, and Yingmin Yi. Inference of transcriptional regulation using multi-objective optimization with gaussian process. In Image and Signal Processing, BioMedical Engineering and Informatics (CISP-BMEI), International Congress on, pages 1820–1825. IEEE, 2016.

## 7 Supplementary Material

### 7.1 Derivatives of a Gaussian process

#### 7.1.1 Definitions

Following Papoulis and Pillai (2002), we define stochastic convergence and stochastic differentiability.

###### Definition

The RV converges to x in the MS sense (limit in mean) if for some x

 limn→∞E(|xn−x|)=0 (20)

###### Definition

The stochastic process x(t) is MS differentiable if for some

 limδt→0E∣∣∣x(t+δt)−x(t)δt−x′(t)∣∣∣=0 (21)

###### Definition

A stochastic process x(t) is called a Gaussian Process, if any finite number of samples of its trajectory are jointly Gaussian distributed according to a previously defined mean function and a covariance matrix, that can be constructed using a predefined kernel function

#### 7.1.2 GP and its derivative are jointly Gaussian

Let .
Let be a Gaussian Process with constant mean and kernel function , assumed to be MS differentiable.

From the definition of GP, we know that and are jointly Gaussian distributed.

 [x(t0)x(t0+δt)]∼N([μμ],Σ) (22)

where using .

Using the linear transformation

 T=1δt[δt0−11] (23)

one can show that

 [x(t0)x(t0+δt)−x(t0)δt]∼N([μ0],TΣTT) (24)

So it is clear that for all , and are jointly Gaussian distributed. Using the assumption that x is differentiable according to the definition in eq. 21 and the fact that convergence in expectation implies convergence in distribution, it is clear that and are jointly Gaussian.

This fact can be easily extended to any finite set of sample times . One can use the exact same procedure to show that the resulting vectors and are jointly Gaussian as well.

#### 7.1.3 Moments of the joint distribution

As shown in the previous section, any finite set of samples is jointly Gaussian together with its derivatives . To calculate the full distribution, it thus suffices to calculate the mean and the covariance between the elements of the full vector

 [x˙x]∼N([μ0],[CϕC′ϕ′CϕC′′ϕ]) (25)

is the predefined kernel matrix of the Gaussian Process.

can be calculated by directly using the linearity of the covariance operator.

 ′Cϕi,j =cov(˙x(ti),x(tj)) (26) =cov(ddax(a)∣∣∣a=ti,x(tj)) =ddacov(x(a),x(tj))a=ti =ddakϕ(a,tj)a=ti

Obviously, is just the transposed of , while can be calculated in exactly the same manner to obtain

#### 7.1.4 Conditional GP for derivatives

To obtain the GP over the derivatives given the states, the joint distribution

 [x˙x]∼N([μ0],[CϕC′ϕ′CϕC′′ϕ]) (28)

has to be transformed. This can be done using standard techniques as described e.g. in section 8.1.3 of Petersen et al. (2008). There, it is written:

Define

 (29)

Then

 p(xb|xa)∼N(^μb,^Σb) (30)

where

 ^μb =μb+ΣTcΣ−1a(xa−μa) (31) ^Σb =Σb−ΣTcΣ−1aΣc (32)

Applied to the above probability distribution, this leads to

 p(˙x|x) ∼N(Dx,A) (33)

using

 D =′CϕCϕ−1 (34) A =C′′ϕ−′CϕCϕ−1C′ϕ (35)

### 7.2 Proof of theorem 1

###### Proof

The proof of this statement follows directly by combining all the previous definitions and marginalizing out all the random variables that are not part of the end result.

First, one starts with the joint density over all variables as stated in equation (15)

 p(x,˙x, y,F1,F2,θ|ϕ,σ,γ)= pGP(x,˙x,y|ϕ,σ)pODE(F1,F2,θ|x,˙x,γ),

where

 pGP(x,˙x,y|ϕ,σ)=p(x|ϕ)p(˙x|x,ϕ)p(y|x,σ)

and

 pODE( F1,F2,θ|x,˙x,γ)= p(θ)p(F1|θ,x)p(F2|˙x,γI)δ(F1−F2).

To simplify this formula, can be reduced by marginalizing out using the properties of the Dirac delta function and the probability density defined in equation (14). The new is then independent of .

 pODE(F1,θ|x,˙x,γ)=p(θ)p(F1|θ,x)N(F1|˙x,γI).

Inserting equation (13) yields

Again, the properties of the Dirac delta function are used to marginalize out . The new is now independent of ,

 pODE(θ|x,˙x,γ)=p(θ)N(f(x,θ)|˙x,γI).

This reduced is now combined with . Observing that the mean and the argument of a normal density are interchangeable and inserting the definition of the GP prior on the derivatives given by equation (5) leads to

 p(x, ˙x,y,θ|ϕ,σ,γ)= p(θ)p(x|ϕ)N(˙x|Dx,A)p(y|x,σ)N(˙x|f(x,θ),γI).

can now be marginalized by observing that the product of two normal densities of the same variable is again a normal density. The formula can be found, e.g., in Petersen et al. (2008). As a result, one obtains

 p(x, y,θ|ϕ,σ,γ)= p(θ)p(x|ϕ)p(y|x,σ)N(f(x,θ)|Dx,A+γI).

It should now be clear that after inserting equations (3) and (4) and renormalizing, we get the final result

 p( x,θ|y,ϕ,γ,σ)∝ p(θ)N(x|0,Cϕ)N(y|x,σ2I)N(f(x,θ)|Dx,A+γI),

concluding the proof of this theorem.

### 7.3 Hyperparameter and kernel selection

As discussed before, the Gaussian process model is defined by a kernel function . For both the hyperparameters and the functional form of there exist many possible choices. Even though the exact choice might not be too important for consistency guarantees in GP regression (Choi and Schervish, 2007), this choice directly influences the amount of observations that are needed for reasonable performance. While there exist some interesting approaches to learn the kernel directly from the data, e.g., Duvenaud et al. (2013) and Gorbach et al. (2017b), these methods can not be applied due to the very low amount of observations of the systems considered in this paper. As in previous approaches, the kernel functional form is thus restricted to simple kernels with few hyperparameters, whose behaviors have already been investigated by the community, e.g., in the kernel cookbook by Duvenaud (2014). Once a reasonable kernel is chosen, it is necessary to fit the hyperparameter and depending on the amount of expert knowledge available, there are different methodologies.

#### 7.3.1 Maximizing the data likelihood

As mentioned e.g. in Rasmussen (2004), it is possible for a Gaussian process model to analytically calculate the marginal likelihood of the observations given the evaluation times and hyperparameters and .

 log( p(y|t,ϕ,σ))= −12yT(Cϕ+σI)−1 −12log|Cϕ+σI| −n2log2π (36)

where is the GPs estimate for the standard deviation of the observation noise and n is the amount of observations.

This equation is completely independent of the ODE system one would like to analyze and depends only on the observations of the states. To fit the GP model to the data, equation (7.3.1) can be maximized w.r.t. and , without incorporating any prior knowledge.

#### 7.3.2 Concurrent optimization

In of Dondelinger et al. (2013), the hyperparameters are not calculated independent of the ODEs. Instead, a prior is defined and their posterior distribution is determined simultaneously with the posterior distribution over states and parameters by sampling from equation (8).

This approach has several drawbacks. As we shall see in section 5, its empirical performance is significantly depending on the hyperpriors. Furthermore, optimizing the joint distribution given equation (8) requires calculating the inverse of the covariance matrices and , which has to be done again and again for each new set of hyperparameters. Due to the computational complexity of matrix inversion, this is significantly slowing down the optimization.

For these reasons, if strong prior knowledge about the hyperparameters is available, it might be better to incorporate it into the likelihood optimization presented in the previous section. There, it could be used as a hyperprior to regularize the likelihood.

#### 7.3.3 Manual tuning

In the variational inference approach Gorbach et al. (2017a), the hyperparameters were assumed to be provided by an expert. If such expert knowledge is available, it should definitely be used since it can improve the accuracy drastically.

### 7.4 Adjusting the GP model

To make more comparable to , the hyperparameters of the kernel must be learned from the data. However, maximizing the data likelihood described in equation (7.3.1) directly using the prior defined in equation (4) will lead to very bad results.

#### 7.4.1 Zero mean observations

The main reason for the poor performance without any pretreatment of is the fact that the zero mean assumption in equation (7.3.1) is a very strong regularization for the amount of data available. As its effect directly depends on the distance of the true values to zero, it will be different for different states in multidimensional systems, further complicating the problem. Thus, it is common in GP regression to manipulate the observations such that they have zero mean.

This procedure can be directly incorporated into the joint density given by equation (8). It should be noted that for multidimensional systems this joint density will factorize over each state k, whose contribution will be given by

 p(xk| yk,θ,ϕ,γ,σ)∝N(~xk|0,Cϕk) ×N(fk(xk,θ)|Dk~xk,Ak+γI) (37)

where

 ~xk =xk−μy,k1

using to denote the mean of the observations of the k-th state and to denote a vector of ones with appropriate length.

It is important to note that this transformation is not just equivalent to exchanging and . While the transformation is not necessary for the observation term, as and would be shifted equally, the original is needed as input to the ODEs. This allows for this transformation without the need to manually account for this in the structural form of the differential equations.

This trick will get rid of some of the bias introduced by the GP prior. In the simulations, this made a difference for all systems, including the most simple one presented in section 5.1.

#### 7.4.2 Standardized states

If the systems get more complex, the states might be magnitudes apart from each other. If one were to use the same hyperparameters for all states, then a deviation would contribute equally to the change in probability, independent of whether the states are of magnitude or . Thus, small relative deviations from the mean of states with large values will lead to stronger changes in the joint probability than large relative deviations of states with small values. This is not a desirable property, which can be partially alleviated by calculating a new set of hyperparameters for each state. However, this problem can be completely nullified by standardizing the data . For the k-th state, this would change its contribution to the joint density to

 p(xk| yk,θ,ϕ,γ,σ)∝N(1σy,k~xk∣∣∣0,Cϕk) (38) ×N(1σy,kyk∣∣∣1σy,kxk,σ2I) ×N(1σy,kfk(xk,θ)∣∣∣1σy,kDk~xk,Ak+γI)

where

 ~xk =xk−μy,k1

and is the standard deviation of the observations of the k-th state.

For complex systems with states on different orders of magnitude, standardization is a must to obtain reasonable performance. Even for the system presented in section 5.2, standardization has a significantly beneficial effect, although the states do not differ greatly in magnitude.

### 7.5 Algorithmic details

For all experiments involving , the R toolbox deGradInfer (Macdonald and Dondelinger, 2017) published alongside Macdonald (2017) was used. For comparability, no priors were used, but the toolbox needed to be supplied with a value for the standard deviation of the observational noise and the true standard deviation of the noise was used. For all other parameters, e.g., amount of chains or samples, the values reported in Dondelinger et al. (2013) were used.

For and , the hyperparameters of the GP model were determined in a preprocessing step identical for both algorithms. After calculating the hyperparameters, the parameters were inferred using the implementation used by Gorbach et al. (2017a).

Both and need to be supplied with , which was treated as a tuning parameter. In principle, this parameter could be found by evaluating multiple candidates in parallel and choosing based on data fit.

For both experiments, the amount of iterations recommended by the toolbox (Macdonald and Dondelinger, 2017) have been used, namely 100’000 iterations for Lotka Volterra and 300’000 iterations for Protein Transduction. This is the same setup that has been used to obtain the parameter estimates shown throughout the paper. The simpler sampling setup of clearly leads to running time savings of about a third compared to . Thus, is not only significantly more accurate, but also significantly faster than . Dondelinger et al. (2013) have shown that this also implies order of magnitude improvements if compared to the running time of approaches based on numerical integration.

All experiments were performed using the ETH cluster.. It should be noted that the algorithms were implemented in different programming languages and the cluster consists of cores with varying computation power, adding some variability to the running time estimates.

### 7.6 Lotka Volterra Figure 9: Example rollout of the Lotka Volterra system showing the state evolution over time. The dots denote the observations while the line represents the ground truth obtained by numerical integration.

As in the previous publications, a squared exponential kernel was used. For and , was set to 0.3, while was provided with the true observation noise standard deviations . The standard deviation of the proposal distribution of was chosen as 0.075 for state proposals and as 0.09 for parameter proposals to roughly achieve an acceptance rate of 0.234. For all algorithms, it was decided to use only one GP to fit both states. This effectively doubles the amount of observations and leads to more stable hyperparameter estimates. As the state dynamics are very similar, this approximation is feasible. Figure 10: Example rollout of the Lotka Volterra system showing the state evolution over time for the high noise case. The dots denote the observations while the line represents the ground truth obtained by numerical integration. (a) AGM (a) AGM (a) AGM

### 7.7 Parameter Distribution Figure 14: Histograms representing the MCMC samples obtained for one example run of the Lotka Volterra system. Each histogram represents the marginal distribution of one ODE parameter.

The MCMC approach of allows to infer the probability distribution over parameters. This is shown for one example rollout in Figure 14. The inferred distributions are close to Gaussian in shape. This likely explains the sampling-like performance of the variational approach , as their assumptions of using a factorized Gaussian proxy distribution over the parameters seems to be a good fit for the true distribution.

### 7.8 Protein Transduction

Figure 15 and Figure 16 show median plots for the states obtained by numerical integration of the inferred parameters of the Protein Transduction system. Figure 15: Median plots for all states of the most difficult benchmark system in the literature, Protein Transduction. The red line is the ground truth, while the black line and the shaded area denote the median and the 75% quantiles of the results of 100 independent noise realizations. FGPGM (middle) is clearly able to find more accurate parameter estimates than AGM (top). Figure 16: Median Plots for all states of the most difficult benchmark system in the literature, Protein Transduction, for the high noise case.The red line is the ground truth, while the black line and the shaded area denote the median and the 75% quantiles of the results of 100 independent noise realizations. FGPGM (middle) is clearly able to find more accurate parameter estimates than AGM (top). Figure 17: Median plots of the numerically integrated states after parameter inference for the FHN system with SNR 10. Ground truth (red), median (black) and 75% quantiles (gray) over 100 independent noise realizations.