In this notebook we consider the latent force model paradigm (Alvarez et al, 2013, Alvarez et al, 2011, Alvarez et al, 2009) to drive spatio-temporal differential equation with a Gaussian process. Latent force models are differential equations whose initial conditions or driving forces are given by a stochastic process. Linear latent force models are linear (partial or ordinary) differential equations. If a linear latent force model is driven by a Gaussian process latent force, this describes a joint Gaussian process distribution across all the variables of interest. The Gaussian process covariance function encodes the relationships between the variables, as proscribed by the differential equation, through the covariance function.
This direction of research is part of our ideas as to how to merge mechanistic and statistical models of data. It also maps onto our ideas about 'Gaussian processes over everything'. The covariance function interelates the protein and mRNA concentrations in the same model.
A model of post-transcriptional processing is formulated to describe the spatio-temporal Drosophila protein expression data (Becker et al, 2013, Alvarez et al, 2011). Protein production is considered to be linearly dependent on the concentration of mRNA at an earlier time point. The model also allows for diffusion of protein between nuclei and linear protein decay. These processes are dependent on the diffusion parameter and the degradation rate of protein respectively.
The coefficients $a$, $b$ and $c$ are unknown. In this study, we use Gaussian process with an exponentiated quadratic kernel as a prior over $y_{x,t}$ (protein). The kernel of $f_{x,t}$ (mRNA) is derived by applying the partial differential operator on the spatial-temporal kernel of protein. The multi-output Gaussian process are developed by combining the covariance matrix of mRNA and protein and their cross covariance.
In [2]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
import pylab as pb
import GPy
import pandas
from pandas import read_csv
The spatio-temporal multi-output partial differential equation covariance functions have been developed with the kernel name GPy.kern.ODE_st in GPy. The inputs are one dimension spatial data, one dimension temporal data and one dimension index which is used to indicate $f$ and $y$.
In [3]:
data = GPy.util.datasets.drosophila_knirps()
The next thing to do is to compose the data ready for presentation to GPy. Here we need to use the time field as the input and a concatanation of the expression levels as the output. We need to provide a corresponding index for to each input to describe what output is being represented.
Now we set up the covariance function with the relevant parameters.
In [4]:
kern = GPy.kern.ODE_st(input_dim=3,
a=1., b=1., c=1.,
variance_Yx=1., variance_Yt=1.,
lengthscale_Yx=15.,
lengthscale_Yt=15.)
With the data correctly presented and the covariance function defined, we are ready to proceed with the Gaussian process regression.
In [5]:
m = GPy.models.GPRegression(data['X'],data['Y'],kern)
The initial value of $a$, $b$ and $c$ are 1. For these choices of covariance function parameters, we can plot the random field of $f$ and $y$ separately.
In [8]:
leng = data['X'].shape[0]
m.plot(fixed_inputs=[(2,0)], which_data_rows = slice(0,leng*2/2))
m.plot(fixed_inputs=[(2,1)], which_data_rows = slice(leng*2/2,leng*2))
Out[8]:
Now we optimize the model, this will take a few minutes.
In [16]:
m.optimize(messages=True)
After optimization, the estimated value of $a$, $b$ and $c$ can be printed.
In [17]:
print m
The plot of the random fields are plotted below.
In [18]:
m.plot(fixed_inputs=[(2,0)], which_data_rows = slice(0,leng*2/2))
pb.savefig("gene.pdf")
m.plot(fixed_inputs=[(2,1)], which_data_rows = slice(leng*2/2,leng*2))
pb.savefig("protein.pdf")
In Kolja Becker's thesis he showed the value of $\alpha$, D and $\lambda$. Should the Gaussian process estimation of $a$, $b$ and $c$ be the same as the value in the thesis? The value of $\lambda$ is equal to 1/$b$.
In [ ]: