Inference of dynamic systems from noisy and sparse data via manifold-constrained Gaussian processes

Following the Bayesian paradigm, we view the $D$-dimensional system $\bm{x}(t)$ to be a realization of the stochastic process $\bm{X}(t)=(X_{1}(t),\ldots,X_{D}(t))$, and the model parameters $\bm{\theta}$ a realization of the random variable $\bm{\Theta}$. In Bayesian statistics, the basis of inference is the posterior distribution, obtained by combining the likelihood function with a chosen prior distribution on the unknown parameters and stochastic processes. Specifically, we impose a general prior distribution $\pi(\cdot)$ on $\bm{\theta}$ and independent GP prior distributions on each component $X_{d}(t)$ so that $X_{d}(t)\sim\mathcal{GP}(\mathcal{\mu}_{d},\mathcal{K}_{d}),~{}t\in[0,T]$, where $\mathcal{K}_{d}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is a positive definite covariance kernel for the GP and $\mathcal{\mu}_{d}:\mathbb{R}\to\mathbb{R}$ is the mean function. Then for any finite set of time points $\bm{\tau}_{d}$, $X_{d}(\bm{\tau}_{d})$ has a multivariate Gaussian distribution with mean vector $\mathcal{\mu}_{d}(\bm{\tau}_{d})$ and covariance matrix $\mathcal{K}_{d}(\bm{\tau}_{d},\bm{\tau}_{d})$. Denote the observations by $\bm{y}({\bm{\tau}})=(\bm{y}_{1}(\bm{\tau}_{1}),\ldots,\bm{y}_{D}(\bm{\tau}_{D}))$, where $\bm{\tau}=(\bm{\tau}_{1},\bm{\tau}_{2},\dots,\bm{\tau}_{D})$ is the collection of all observation time points and each component $X_{d}$ can have its own set of observation times $\bm{\tau}_{d}=(\tau_{d,1},\ldots,\tau_{d,N_{d}})$. If the $d$-th component is not observed, then $N_{d}=0$, and $\bm{\tau}_{d}=\emptyset$. $N=N_{1}+\cdots+N_{D}$ is the total number of observations. We note that for the remainder of the paper, the notation $t$ shall refer to time generically, while $\tau$ shall refer specifically to the observation time points.

As an illustrative example, consider the dynamic system in [1] that governs the oscillation of Hes1 mRNA ($M$) and Hes1 protein ($P$) levels in cultured cells, where it is postulated that a Hes1-interacting ($H$) factor contributes to a stable oscillation, a manifestation of biological rhythm [2]. The ODEs of the three-component system $X=(P,M,H)$ are

$\displaystyle\mathbf{f}(X,\bm{\theta},t)=\begin{pmatrix}-aPH+bM-cP\\ -dM+\frac{e}{1+P^{2}}\\ -aPH+\frac{f}{1+P^{2}}-gH\end{pmatrix},$

where $\bm{\theta}=(a,b,c,d,e,f,g)$ are the associated parameters. In Fig 1 (left panel) we show noise-contaminated data generated from the system, which closely mimics the experimental setup described in [1]: $P$ and $M$ are observed at 15-minute intervals for 4 hours but $H$ is never observed. In addition, $P$ and $M$ observations are asynchronous: starting at time 0, every 15 minutes we observe $P$; starting at 7.5 minutes, every 15 minutes we observe $M$; $P$ and $M$ are never observed at the same time. It can be seen that the mRNA and protein levels exhibit the behavior of regulation via negative feedback.

The goal here is to infer the seven parameters of the system: $a,b$ govern the rate of protein synthesis in the presence of the interacting factor; $c,d,g$ are the rates of decomposition; and $e,f$ are inhibition rates. The unobserved $H$ component poses a challenge for most existing methods that do not use numerical integration, but is capably handled by MAGI: the $P$ and $M$ panels of Fig 1 show that our inferred trajectories provide good fits to the observed data, and the $H$ panel shows that the dynamics of the entirely unobserved $H$ component are largely recovered as well. We emphasize that these trajectories are inferred without any use of numerical solvers. We shall return to the Hes1 example in detail in the Results section.

Intuitively, the GP prior on $\bm{X}(t)$ facilitates computation as GP provides closed analytical forms for $\bm{\dot{X}}(t)$ and $\bm{X}(t)$, which could bypass the need for numerical integration. In particular, with a GP prior on $\bm{X}(t)$, the conditional distribution of $\bm{\dot{X}}(t)$ given $\bm{X}(t)$ is also a GP with its mean function and covariance kernel completely specified. This GP specification for the derivatives $\bm{\dot{x}}(t)$, however, is inherently incompatible with the ODE model because (1) also completely specifies $\bm{\dot{x}}(t)$ given $\bm{x}(t)$ (via the function $\mathbf{f}$). As a key novel contribution of our method, MAGI addresses this conceptual incompatibility by constraining the GP to satisfy the ODE model in (1). To do so, we first define a random variable $W$ quantifying the difference between stochastic process $\bm{X}(t)$ and the ODE structure with a given value of the parameter $\bm{\theta}$:

$$W=\sup_{t\in[0,T],d\in\{1,\ldots,D\}}|\dot{X}_{d}(t)-\mathbf{f}(\bm{X}(t),\bm{\theta},t)_{d}|.$$

(2)

$W\equiv 0$ if and only if ODEs with parameter $\bm{\theta}$ are satisfied by $\bm{X}(t)$. Therefore, ideally the posterior distribution for $\bm{X}(t)$ and $\bm{\theta}$ given the observations $\bm{y}(\bm{\tau})$ and the ODE constraint, $W\equiv 0$, is (informally)

$$p_{\bm{\Theta},\bm{X}(t)|W,\bm{Y(\tau})}(\bm{\theta},\bm{x}(t)|W=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau})).$$

(3)

While (3) is the ideal posterior, in reality $W$ is not generally computable. In practice we approximate $W$ by finite discretization on the set $\bm{I}=(t_{1},t_{2},\ldots,t_{n})$ such that $\bm{\tau}\subset\bm{I}\subset[0,T]$ and similarly define $W_{\bm{I}}$ as

$$W_{\bm{I}}=\max_{t\in\bm{I},d\in\{1,\ldots,D\}}|\dot{X}_{d}(t)-\mathbf{f}(\bm{X}(t),\bm{\theta},t)_{d}|.$$

(4)

Note that $W_{\bm{I}}$ is the maximum of a finite set, and $W_{\bm{I}}\to W$ monotonically as $\bm{I}$ becomes dense in $[0,T]$. Therefore, the practically computable posterior distribution is

$$p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I},\bm{Y(\tau})}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau})),$$

which is the joint conditional distribution of $\bm{\theta}$ and $\bm{X}(\bm{I})$ together. Thus, effectively, we simultaneously infer both the parameters and the unobserved trajectory $\bm{X}(\bm{I})$ from the noisy observations $\bm{y}(\bm{\tau})$.

Under Bayes’ rule, we have

	$\displaystyle p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I},\bm{Y(\tau})}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau}))$
	$\displaystyle\propto P({\bm{\Theta}}=\bm{\theta},\bm{X}(\bm{I})=\bm{x}(\bm{I}),W_{\bm{I}}=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau})),$

where the right hand side can be decomposed as

	$\displaystyle P({\bm{\Theta}}=\bm{\theta},\bm{X}(\bm{I})=\bm{x}(\bm{I}),W_{\bm{I}}=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau}))$
	$\displaystyle=\pi_{\bm{\Theta}}(\bm{\theta})\times\underbrace{P(\bm{X}(\bm{I})=\bm{x({I})}|{\bm{\Theta}}=\bm{\theta})}_{(1)}$
	$\displaystyle\qquad\times\underbrace{P(\bm{Y(\bm{\tau})}=\bm{y(\bm{\tau})}|\bm{X}(\bm{I})=\bm{x}(\bm{I}),{\bm{\Theta}}=\bm{\theta})}_{(2)}$
	$\displaystyle\qquad\times\underbrace{P(W_{\bm{I}}=0|\bm{Y(\bm{\tau})}=\bm{y(\bm{\tau})},\bm{X}(\bm{I})=\bm{x}(\bm{I}),{\bm{\Theta}}=\bm{\theta})}_{(3)}.$

The first term (1) can be simplified as $P(\bm{X}(\bm{I})=\bm{x({I})}|{\bm{\Theta}}=\bm{\theta})=P(\bm{X}(\bm{I})=\bm{x({I})})$ due to the prior independence of $\bm{X}(\bm{I})$ and $\bm{\Theta}$; it corresponds to the GP prior on $\bm{X}$. The second term (2) corresponds to the noisy observations. The third term (3) can be simplified as

	$\displaystyle P(W_{\bm{I}}=0|\bm{Y(\bm{\tau})}=\bm{y(\bm{\tau})},\bm{X}(\bm{I})=\bm{x}(\bm{I}),{\bm{\Theta}}=\bm{\theta})$
	$\displaystyle=P(\bm{\dot{X}({I})}-\mathbf{f}(\bm{x}(\bm{I}),\bm{\theta},t_{\bm{I}})=\bm{0}|\bm{Y(\bm{\tau})}=\bm{y(\bm{\tau})},\bm{X}(\bm{I})=\bm{x}(\bm{I}),{\bm{\Theta}}=\bm{\theta})$
	$\displaystyle=P(\bm{\dot{X}({I})}-\mathbf{f}(\bm{x}(\bm{I}),\bm{\theta},t_{\bm{I}})=\bm{0}|\bm{X}(\bm{I})=\bm{x({I})})$
	$\displaystyle=P(\bm{\dot{X}({I})}=\mathbf{f}(\bm{x}(\bm{I}),\bm{\theta},t_{\bm{I}})|\bm{X}(\bm{I})=\bm{x({I})}),$

which is the conditional density of $\bm{\dot{X}({I})}$ given $\bm{X}(\bm{I})$ evaluated at $\mathbf{f}(\bm{x}(\bm{I}),\bm{\theta},t_{\bm{I}})$. All three terms are multivariate Gaussian: the third term is Gaussian because $\bm{\dot{X}({I})}$ given $\bm{X}(\bm{I})$ has a multivariate Gaussian distribution as long as the kernel $\mathcal{K}$ is twice differentiable.

Therefore, the practically computable posterior distribution simplifies to

	$\displaystyle p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I},\bm{Y(\tau})}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0,\bm{Y}(\bm{\tau})=\bm{y}(\bm{\tau}))$		(5)
	$\displaystyle\propto\pi_{\bm{\Theta}}(\bm{\theta})\exp\Big{\{}-\frac{1}{2}\sum_{d=1}^{D}\Big{[}$
	$\displaystyle+\underbrace{|\bm{I}|\log(2\pi)+\log{|C_{d}|}+\left\|x_{d}(\bm{I})-\mu_{d}(\bm{I})\right\|_{C_{d}^{-1}}^{2}}_{(1)}$
	$\displaystyle+\underbrace{|\bm{I}|\log(2\pi)+\log{|K_{d}|}+\left\|\mathbf{f}_{d,\bm{I}}^{\bm{x},\bm{\theta}}-\dot{\mu}_{d}(\bm{I})-m_{d}\{x_{d}(\bm{I})-\mu_{d}(\bm{I})\}\right\|_{K_{d}^{-1}}^{2}}_{(3)}$
	$\displaystyle+\underbrace{N_{d}\log(2\pi\sigma_{d}^{2})+\left\|x_{d}(\bm{\tau}_{d})-y_{d}(\bm{\tau}_{d})\right\|_{\sigma_{d}^{-2}}^{2}}_{(2)}\Big{]}\Big{\}}$

where $\|\bm{v}\|_{A}^{2}=\bm{v}^{\intercal}A\bm{v}$, $|\bm{I}|$ is the cardinality of $\bm{I}$, $\mathbf{f}_{d,\bm{I}}^{\bm{x},\bm{\theta}}$ is short for the $d$-th component of $\mathbf{f}(\bm{x}(\bm{I}),\bm{\theta},t_{\bm{I}})$, and the multivariate Gaussian covariance matrix $C_{d}$ and the matrix $K_{d}$ can be derived as follows for each component $d$:

$$\begin{cases}C&=\mathcal{K}(\bm{{I}},\bm{I})\\ m&=\mathcal{{}^{\prime}K}(\bm{I},\bm{I})\mathcal{K}(\bm{I},\bm{I})^{-1}\\ K&=\mathcal{K^{\prime\prime}}(\bm{I},\bm{I})-\mathcal{{}^{\prime}K}(\bm{I},\bm{I})\mathcal{K}(\bm{I},\bm{I})^{-1}\mathcal{K^{\prime}}(\bm{I},\bm{I})\end{cases}$$

(6)

where $\mathcal{{}^{\prime}K}=\frac{\partial}{\partial s}\mathcal{K}(s,t)$, $\mathcal{K^{\prime}}=\frac{\partial}{\partial t}\mathcal{K}(s,t)$, and $\mathcal{K^{\prime\prime}}=\frac{\partial^{2}}{\partial s\partial t}\mathcal{K}(s,t)$.

In practice we choose the Matern kernel $\mathcal{K}(s,t)=\phi_{1}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{l}{\phi_{2}}\right)^{\nu}B_{\nu}\left(\sqrt{2\nu}\frac{l}{\phi_{2}}\right)$ where $l=|s-t|$, $\Gamma$ is the Gamma function and $B_{\nu}$ is the modified Bessel function of the second kind, and the degree of freedom $\nu$ is set to be 2.01 to ensure that the kernel is twice differentiable. $\mathcal{K}$ has two hyper-parameters $\phi_{1}$ and $\phi_{2}$. Their meaning and specification are discussed in the Materials and Methods section.

With the posterior distribution specified in (5), we use Hamiltonian Monte Carlo (HMC) [7] to obtain samples of $\bm{X_{I}}$ and the parameters together. At the completion of HMC sampling, we take the posterior mean of $\bm{X_{I}}$ as the inferred trajectory, and the posterior means of the sampled parameters as the parameter estimates. Throughout the MAGI computation, no numerical integration is ever needed.

The problem of dynamic system inference has been studied in the literature, which we now briefly review. We first note that a simple approach to constructing the ‘ideal’ likelihood function is according to $p(\bm{y}(\bm{t})|\hat{\bm{x}}(\bm{t},\bm{\theta},\bm{x}(0)),\sigma)$, where $\hat{\bm{x}}(\bm{t},\bm{\theta},\bm{x}(0))$ is the numerical solution of the ODE obtained by numerical integration given $\bm{\theta}$ and the initial conditions. This approach suffers from a high computational burden: numerical integration is required for every $\bm{\theta}$ sampled in an optimization or Markov chain Monte Carlo (MCMC) routine [8]. Smoothing methods have been useful for eliminating the dependence on numerical ODE solutions, and an innovative penalized likelihood approach [9] uses a B-spline basis for constructing estimated functions to simultaneously satisfy the ODE system and fit the observed data. While in principle the method in [9] can handle an unobserved system component, substantive manual input is required as we show in the Results, which contrasts with the ready-made solution that MAGI provides.

As an alternative to the penalized likelihood approach, GPs are a natural candidate for fulfilling the smoothing role in a Bayesian paradigm due to their flexibility and analytic tractability [10]. The use of GPs to approximate the dynamic system and facilitate computation has been previously studied by a number of authors [8, 11, 12, 13, 14, 15]. The basic idea is to specify a joint GP over $\bm{y},\bm{x},\dot{\bm{x}}$ with hyperparameters $\phi$, and then provide a factorization of the joint density $p(\bm{y},\bm{x},\dot{\bm{x}},\bm{\theta},\phi,\sigma)$ that is suitable for inference. The main challenge is to find a coherent way to combine information from two distinct sources: the approximation to the system by the GP governed by hyperparameters $\phi$, and the actual dynamic system equations governed by parameters $\bm{\theta}$. In [8, 11], the factorization proposed is $p(\bm{y},\bm{x},\dot{\bm{x}},\bm{\theta},\phi,\sigma)=p(\bm{y}|\bm{x},\sigma)p(\dot{\bm{x}}|\bm{x},\bm{\theta},\phi)p(\bm{x}|\phi)p(\phi)p(\bm{\theta})$, where $p(\bm{y}|\bm{x},\sigma)$ comes from the observation model and $p(\bm{x}|\phi)$ comes from the GP prior as in our approach. However, there are significant conceptual difficulties in specifying $p(\dot{\bm{x}}|\bm{x},\bm{\theta},\phi)$: on one hand, the distribution of $\dot{\bm{x}}$ is completely determined by the GP given $\bm{x}$, while on the other hand $\dot{\bm{x}}$ is completely specified by the ODE system $\dot{\bm{x}}=\mathbf{f}(\bm{x},\bm{\theta},t)$; these two are incompatible. Previous authors have attempted to circumvent this incompatibility of the GP and ODE system: [8, 11] use a product-of-experts heuristic by letting $p(\dot{\bm{x}}|\bm{x},\bm{\theta},\phi)\propto p(\dot{\bm{x}}|\bm{x},\phi)p(\dot{\bm{x}}|\bm{x},\bm{\theta})$, where the two distributions in the product come from the GP and a noisy version of the ODE, respectively. In [15], the authors arrive at the same posterior as [8, 11] by assuming an alternative graphical model that bypasses the product of experts heuristic; nonetheless, the method requires working with an artificial noisy version of the ODE. In [12], the authors start with a different factorization: $p(\bm{y},\bm{x},\dot{\bm{x}},\bm{\theta},\phi,\sigma)=p(\bm{y}|\dot{\bm{x}},\phi,\sigma)p(\dot{\bm{x}}|\bm{x},\bm{\theta})p(\bm{x}|\phi)p(\phi)p(\bm{\theta})$, where $p(\bm{y}|\dot{\bm{x}},\phi)$ and $p(\bm{x}|\phi)$ are given by the GP and $p(\dot{\bm{x}}|\bm{x},\bm{\theta})$ is a Dirac delta distribution given by the ODE. However, this factorization is incompatible with the observation model $p(\bm{y}|\bm{x},\sigma)$ as discussed in detail in [16]. There is other related work that uses GPs in an ad hoc partial fashion to aid inference. In [13], GP regression is used to obtain the means of $\bm{x}$ and $\dot{\bm{x}}$ for embedding within an Approximate Bayesian Computation estimation procedure. In [14], GP smoothing is used during an initial burn-in phase as a proxy for the likelihood, before switching to the ‘ideal’ likelihood to obtain final MCMC samples. While empirical results from the aforementioned studies are promising, a principled statistical framework for inference that addresses the previously noted conceptual incompatibility between the GP and ODE specifications is lacking. Our work presents one such principled statistical framework through the explicit manifold constraint. MAGI is therefore distinct from recent GP-based approaches [11, 15] or any other Bayesian analogs of [9].

In addition to the conceptual incompatibility, none of the existing methods that do not use numerical integration offer a practical solution for a system with unobserved component(s), which highlights another unique and important contribution of our approach.

The Hes1 model described in the Introduction demonstrates inference on a system with an unobserved component and asynchronous observation times. This section continues the inference of this model. Ref [1] studied the theoretical oscillation behavior using parameter values $a=0.022$, $b=0.3$, $c=0.031$, $d=0.028$; $e=0.5$, $f=20$, $g=0.3$, which leads to one oscillation cycle approximately every 2 hours. Ref [1] also set the initial condition at the lowest value of $P$ when the system is in oscillation equilibrium [1]: $P=1.439$, $M=2.037$, $H=17.904$. The noise level in our simulation is derived from [1] where the standard error based on repeated measures are reported to be around 15% of the $P$ (protein) level and $M$ (mRNA) level, so we set the simulation noise to be multiplicative following a log-normal distribution with standard deviation 0.15, and throughout this example we assume the noise level $\sigma$ is known to be 0.15 from repeated measures reported in [1]. The $H$ component is never observed. Owing to the multiplicative error on the strictly positive system, we apply our method to the log-transformed ODEs, so that the resulting error distributions are Gaussian. To the best of our knowledge, MAGI is the only one that provides a practical and complete solution for handling unobserved component cases like this example.

We generate 2000 simulated datasets based on the above setup for the Hes1 system. The left-most panel in Fig 1 shows one example dataset. For each dataset, we use MAGI to infer the trajectories and estimate the parameters. We use the posterior mean of $X_{t}=(P,M,H)_{t}$ as the inferred trajectories for the system components, which are generated by MAGI without using any numerical solver. Fig 1 summarizes the inferred trajectories across the 2000 simulated datasets, showing the median of the inferred trajectories of $X_{t}$ together with the 95% interval of the inferred trajectories represented by the 2.5% and 97.5% percentiles. The posterior mean of $\bm{\theta}=(a,b,c,d,f,e,g)$ is our estimate of the parameters. Table 1 summarizes the parameter estimates across the 2000 simulated datasets, by showing their means and standard deviations. Fig 1 shows that MAGI recovers the system well, including the completely unobserved $H$ component. Table 1 shows that MAGI also recovers the system parameters well, except for the parameters that only appear in the equation for the unobserved $H$ component, which we will discuss shortly. Together, Fig 1 and Table 1 demonstrate that MAGI can recover the entire system without any usage of a numerical solver, even in the presence of unobserved component(s).

To further assess MAGI, we compare with two methods: Adaptive Gradient Matching (AGM) of Ref [11] and Fast Gaussian process based Gradient Matching (FGPGM) of Ref [15], representing the state-of-the-art of inference methods based on GPs. For fair comparison, we use the same benchmark systems, scripts and software provided by the authors for performance assessment, and run the software using the settings recommended by the authors. The benchmark systems include the FitzHugh-Nagumo (FN) equations [17] and a protein transduction model [18].

We begin by summarizing the computational scheme of MAGI. Overall, we use Hamiltonian Monte Carlo (HMC) [7] to obtain samples of $\bm{X_{I}}$ and the parameters from their joint posterior distribution. Details of the HMC sampling are included in the SI section ‘Hamiltonian Monte Carlo’. At each iteration of HMC, $\bm{X_{I}}$ and the parameters²²2The parameters here refer to $\bm{\theta}$ and $\sigma$. If the noise level $\sigma$ is known a priori, the parameters then refer to $\bm{\theta}$ only. are updated together with a joint gradient, with leapfrog step sizes automatically tuned during the burn-in period to achieve an acceptance rate between 60-90%. At the completion of HMC sampling (and after discarding an appropriate burn-in period for convergence), we take the posterior means of $\bm{X_{I}}$ as the inferred trajectories, and the posterior means of the sampled parameters as the parameter estimates. The techniques we use to temper the posterior and speed up the computations are discussed in the following ‘Prior tempering’ subsection and ‘Techniques for computational efficiency’ in the SI.

Several steps are taken to initialize the HMC sampler. First, we apply a GP fitting procedure to obtain values of $\bm{\phi}$ and $\sigma$ for the observed components; the computed values of the hyper-parameters $\bm{\phi}$ are subsequently held fixed during the HMC sampling, while the computed value of $\sigma$ is used as the starting value in the HMC sampler. (If $\sigma$ is known, the GP fitting procedure is used to obtain values of $\bm{\phi}$ only.) Second, starting values of $\bm{X_{I}}$ for the observed components are obtained by linearly interpolating between the observation time points. Third, starting values for the remaining quantities – $\bm{\theta}$ and $(\bm{X_{I}},\bm{\phi})$ for any unobserved component(s) – are obtained by optimization of the posterior as described below.

The GP prior $X_{d}(t)\sim\mathcal{GP}(\mathcal{\mu}_{d},\mathcal{K}_{d}),~{}t\in[0,T]$, is on each component $X_{d}(t)$ separately. The Gaussian process Matern kernel $\mathcal{K}(l)=\phi_{1}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{l}{\phi_{2}}\right)^{\nu}B_{\nu}\left(\sqrt{2\nu}\frac{l}{\phi_{2}}\right)$ has two hyper-parameters that are held fixed during sampling: $\phi_{1}$ controls overall variance level of the GP, while $\phi_{2}$ controls the bandwidth for how much neighboring points of the GP affect each other.

When the observation noise level $\sigma$ is unknown, values of $(\phi_{1},\phi_{2},\sigma)$ are obtained jointly by maximizing GP fitting without conditioning on any ODE information, namely:

	$\displaystyle(\tilde{\bm{\phi}},\tilde{\sigma})$	$\displaystyle=\operatorname*{arg\,max}_{\bm{\phi},\sigma}p(\bm{\phi},\sigma^{2}|\bm{y_{{}_{I_{0}}}})$
		$\displaystyle=\operatorname*{arg\,max}_{\bm{\phi},\sigma}\pi_{\Phi_{1}}(\phi_{1})\pi_{\Phi_{2}}(\phi_{2})\pi_{\sigma}(\sigma^{2})p(\bm{y_{{}_{I_{0}}}}|\bm{\phi},\sigma^{2})$		(7)

where $\bm{y_{{}_{I_{0}}}}|\bm{\phi},\sigma\sim\mathcal{N}(0,\mathcal{K}_{\phi}+\sigma^{2})$. The index set $I_{0}$ is the smallest evenly spaced set such that all observation time points in this component are in $I_{0}$, i.e., $\bm{\tau}\subseteq I_{0}$. The priors $\pi_{\Phi_{1}}(\phi_{1})$ and $\pi_{\sigma}(\sigma^{2})$ for the variance parameter $\phi_{1}$ and $\sigma$ are set to be flat. The prior $\pi_{\Phi_{2}}(\phi_{2})$ for the bandwidth parameter $\phi_{2}$ is set to be a Gaussian distribution: (a) the mean $\mu_{\Phi_{2}}$ is set to be half of the period corresponding to the frequency that is the weighted average of all the frequencies in the Fourier transform of $y$ on $I_{0}$ (the values of $y$ on $I_{0}$ are linearly interpolated from the observations at $\bm{\tau}$), where the weight on a given frequency is the squared modulus of the Fourier transform with that frequency, and (b) the standard deviation is set such that $T$ is three standard deviations away from $\mu_{\Phi_{2}}$. This Gaussian prior on $\phi_{2}$ serves to prevent it from being too extreme. In the subsequent sampling of $(\bm{\theta},\bm{X_{\tau}},\sigma^{2})$, the hyper-parameters $\bm{\phi}$ are fixed at $\tilde{\bm{\phi}}$ while $\tilde{\sigma}$ gives the starting value of $\sigma$ in the HMC sampler.

If $\sigma$ is known, then values of $(\phi_{1},\phi_{2})$ are obtained by maximizing

$$\tilde{\bm{\phi}}=\operatorname*{arg\,max}_{\bm{\phi}}p(\bm{\phi}|\bm{y_{{}_{I_{0}}}},\sigma^{2})=\operatorname*{arg\,max}_{\bm{\phi}}\pi_{\Phi_{1}}(\phi_{1})\pi_{\Phi_{2}}(\phi_{2})p(\bm{y_{{}_{I_{0}}}}|\bm{\phi},\sigma^{2})$$

(8)

and held fixed at $\tilde{\bm{\phi}}$ in the subsequent HMC sampling of $(\bm{\theta},\bm{X_{\tau}})$. The priors for $\phi_{1}$ and $\phi_{2}$ are the same as previously defined.

To provide starting values of $\bm{X_{I}}$ for the HMC sampler, we use the values of $\bm{Y_{\tau}}$ at the observation time points and linearly interpolate the remaining points in $\bm{I}$.

To provide starting values of $\bm{\theta}$ for the HMC sampler, we optimize the posterior (5) as a function of $\bm{\theta}$ alone, holding $\bm{X_{I}}$ and $\sigma$ unchanged at their starting values, when there is no unobserved component(s). The optimized $\bm{\theta}$ is then used as the starting value for the HMC sampler in this case.

Separate treatment is needed for the setting of $\bm{\phi}$ and initialization of $(\bm{\theta}$, $\bm{X_{I}})$ for the unobserved component(s). We use an optimization procedure that seeks to maximize the full posterior in (5) as a function of $\bm{\theta}$ together with $\bm{\phi}$ and the whole curve of $\bm{X_{I}}$ for unobserved components, while holding the $\sigma$, $\bm{\phi}$ and $\bm{X_{I}}$ for the observed components unchanged at their initial value discussed above. We thereby set $\bm{\phi}$ for the unobserved component, and the starting values of $\bm{\theta}$ and $\bm{X_{I}}$ for unobserved components at the optimized value. In the subsequent sampling, the hyper-parameters are fixed at the optimized $\bm{\phi}$, while the HMC sampling starts at the $\bm{\theta}$ and the $\bm{X_{I}}$ obtained by this optimization.

After $\bm{\phi}$ is set, we use a tempering scheme to control the influence of the GP prior relative to the likelihood during HMC sampling. Note that (5) can be written as

$\displaystyle\begin{split}&p_{{\bm{\Theta}},\bm{X}(\bm{I})|\bm{Y(\tau)},W_{I}}(\bm{\theta},\bm{x}(\bm{I})|\bm{y}(\bm{\tau}),W_{\bm{I}}=0)\\ \propto&p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I}}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0)p_{\bm{Y(\tau)}|\bm{X}(\bm{\tau})}(\bm{y}(\bm{\tau})|\bm{x}(\bm{\tau})).\end{split}$

(9)

As the cardinality of $|\bm{I}|$ increases with more discretization points, the prior part $p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I}}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0)$ grows, while the likelihood part $p_{\bm{Y(\tau)}|\bm{X}(\bm{\tau})}(\bm{y}(\bm{\tau})|\bm{x}(\bm{\tau}))$ stays unchanged. Thus, to balance the influence of the prior, we introduce a tempering hyper-parameter $\beta$ with the corresponding posterior

$\displaystyle\begin{split}&p_{{\bm{\Theta}},\bm{X_{I}}|W_{I},\bm{Y_{\tau}}}^{(\beta)}(\bm{\theta},\bm{x_{I}}|0,\bm{y_{\tau}})\\ \propto&p_{{\bm{\Theta}},\bm{X}(\bm{I})|W_{I}}(\bm{\theta},\bm{x}(\bm{I})|W_{\bm{I}}=0)^{1/\beta}p_{\bm{Y(\tau)}|\bm{X}(\bm{I})}(\bm{y}(\bm{\tau})|\bm{x}(\bm{I}))\\ \propto&\pi_{\bm{\Theta}}(\bm{\theta})\exp\Big{\{}-\frac{1}{2}\sum_{d=1}^{D}\Big{[}N_{d}\log(2\pi\sigma_{d}^{2})+\left\|(x_{d}(\bm{\tau}_{d})-y_{d}(\bm{\tau}_{d}))\right\|_{\sigma_{d}^{-2}}^{2}\\ &\qquad+\frac{1}{\beta}\Big{(}\left\|x_{d}(\bm{I})-\mu_{d}(\bm{I})\right\|_{C_{d}^{-1}}^{2}\\ &\qquad\qquad\quad+\left\|\mathbf{f}_{d,\bm{I}}^{\bm{x},\bm{\theta}}-\dot{\mu}_{d}(\bm{I})-m_{d}(x_{d}(\bm{I})-\mu_{d}(\bm{I}))\right\|_{K_{d}^{-1}}^{2}\Big{)}\Big{]}\Big{\}}\end{split}$

A useful setting that we recommend is $\beta={D|\bm{I}|}/N$, where $D$ is the number of system components, $|\bm{I}|$ is the number of discretization time points, and $N=\sum_{d=1}^{D}N_{d}$ is the total number of observations. This setting aims to balance the likelihood contribution from the observations with the total number of discretization points.

All of the data used in the article are simulation data. The details, including the models to generate the simulation data, are described in Results and the SI. Our software package also includes complete replication scripts for all the data and examples.

The research of S.W.K.W. is supported in part by Discovery Grant RGPIN-2019-04771 from the Natural Sciences and Engineering Research Council of Canada. The research of S.C.K. is supported in part by NSF Grant DMS-1810914.