Approximate Inference for partly observed
Continuous Time Markov Processes
Manfred Opper

Computer Sciencecollaboration with:Andreas Ruttor, Florian Stimberg (TU Berlin)
Dan Cornford, Yuan Shen, Remi Barillec, Michail Vrettas (Aston U)
John Shawe–Taylor, Cédric Archambeau (UCL)
Guido Sanguinetti (Edinburgh U)

Overview• Examples of partly observed continuous time Markov processes• Reminder of Markov processes• Monte Carlo, exact and linear noise approximations• Variational approaches to inference• Applications to a model of transcriptional regulation• OutlookStochastic Lotka Volterra Model2 interacting species (animals): predators and prey.
Numbers of animals: XPrey and XPred.

XPrey → XPrey + 1 with Rate αXPrey
XPrey → XPrey − 1 with Rate βXPreyXPred
XPred → XPred + 1 with Rate δXPreyXPred
XPred → XPred − 1 with Rate γXPredA realisation of the stochastic process and the reaction constants

Discrete noisy observations from the continuous time series

What we would like to predict ...

Gene expression

Simple model of autoregulatory network
2 interacting molecules: mRNA and a Protein
Number of mRNA and Protein molecules: XmRNA, XProt

X = (XmRNA, XProt)

XProt
XProt
XmRNA
XmRNA→→→→XProt + 1
XProt − 1
XmRNA − 1
XmRNA + 1::::

with Ratewith Ratewith Ratewith RateγXmRNA ,δXProt ,βXmRNA ,α(1 − αcΘ(XProt − θc)) ,where Θ(x) = 1 if x ≥ 0 and Θ(x) = 0 for x < 0.

Simulation of process

(Noisy) observations at discrete times.

Inference of state variable.

Motion in double-well potential

dX = f (X)dt + σ 2dW.
(x)with f (x) = − dVdx

and V (x) is a double well potential

A sample path might look like thisObservations & optimal prediction

If we forget about the model and use a Bayesian
regression ...

A reminder of Markov Processes• Discrete time (Markov chain)
P (Xn+1 ∈ A|Xn = x, Xn−1 = xn−1, . . . , X0 = x0)= P (Xn+1 ∈ A|Xn = x)

• Continuous time
P (Xt+dt ∈ A|X0:t = x0:t) = P (Xt+dt ∈ A|Xt = xt)
where x0:t = {x(s)|0 ≤ s ≤ t} denotes the entire path.

Jump Processes
Assume that Xt jumps between discrete states.

• Short time behaviour of transition kernel defined by transition ratef:
P (y, t + ∆t|x, t) ' δxy + f (y|x, t)∆t
Pfor ∆t → 0. Note that f (x|x, t) = − z6=x f (z|x, t).

• Density for time t to next event (homogeneous process)p(t) = −f (x|x) ef (x,x)t

Discrete event simulation
1. initialise X0.
2. Let Xt = x. Simulate waiting time t0 to next event with density0p(t0) = −f (x|x) ef (x|x)t .
3. Put t := t + t00f (x |x)
4. Draw new state x0 with prob − f (x|x)

5. Output Xt = x0

6. If t < Tmax return to step 2.

This is the basis for the Gillespie algorithm.

ExamplesPoisson Process: Xt ∈ N0 and f (x0|x) = λδx0,x+1

20x(t)

15

10

500

2

46t

810

Telegraph Process: Xt ∈ {0, 1}
f (x0|x) = λδx0,1−xx(t)

100

2

46t

810

Forward and Backward equations• Forward (Master) equation iX h dP (x, t|x0, t0)000=P x , t|x0, t0 f x|x − P (x, t|x0, t0) f x |x .
dt
x0 6=x• Kolmogorov Backward equationnoX
dP (x0, t0|x, t)
0000=f (y|x) P (x , t |x, t) − P (x , t |y, t)
dty6=xDiffusion approximation.
Assume that jumps ∆x = x0 − x << x. Try approximation by a continuous state Markov process
Introduce drift (1. Jump moment)i
1 hE Xt+∆t − Xt|Xt = x =f (x) = lim∆tX
f (x0|x)(x0 − x) ,x0 6=x

and diffusion matrix (2. jump moment)i1 h>D(x) = lim
E (Xt+∆t − Xt)(Xt+∆t − Xt) |Xt = x =∆tX(x0 − x)f (x0|x)(x0 − x)> .
x0 6=xDrift & diffusion for Lotka - Volterrah(X1 → X1 + 1) = αX1 ,h(X1 → X1 − 1) = βX1X2 ,h(X2 → X2 + 1) = δX1X2 , h(X2 → X2 − 1) = γX2 ,Deterministic limit
Neglect fluctuations (diffusion), i.e. assume
Xt+∆t − Xt = f (Xt)∆t.
Hence, Xt ≈ m(t) = E[Xt]. This yields ’classical’ rate equationdm≈ f (m(t))
dtDiffusion approximation
Xt+∆t − Xt = f (Xt)∆t + D1/2(Xt)√∆t k .

with k i.i.d. ∼ N (0, I), has the same 2 conditional moments f (x) andD(x). Jump sizes will go → 0 and process becomes continuous!
We will write this (in the limit) as an (Ito) stochastic differentialequation
dXt = f (Xt)dt + D1/2(Xt)dWt

Short time behaviour of transition kernelfor Ito SDE

0p x , t + ∆t|x, t =>
1  01−10exp −(x − x − f (x)∆t D (x) x − x − f (x)∆t=|2πD(x)|2∆tSimulations of Jump and corresponding diffusion
processes
protein300

molecules

250

200

150

100
0

0,51

1,5t

2

2,53

Forward and backward equations for diffusion
processes.
Let x = (x1, . . . xd) and f (x) = (f1(x), . . . , fd(x))• Forward (Fokker - Planck) equation∂p(x, t|x0, t0)=∂t
1 X ∂2 −(fi(x)p(x, t|x0, t0)) +Dij (x)p(x, t|x0, t0)∂x2∂x∂xii jiijX ∂• Kolmogorov Backward equation2XX∂p(x0, t0|x, t)∂1∂=−fi(x)p(x0, t0|x, t)−
Dij (x)p(x0, t0|x, t)
∂t∂xi
2 ij
∂xi∂xji

Short summary• Markov jump processes: Discrete states X00
p X , t + ∆t|X, t ' δX 0X + ∆t f X |X
for ∆t → 0• Diffusion processes: (Ito) stochastic differential equation for stateX(t) ∈ Rd
1/2 (X(t)) ×dX(t) = f(X(t))
dt+DdW(t)| {z }|{z}| {z }
DriftDiffusion
Wiener process

Limit of discrete time process Xk∆Xk ≡ Xk+1 − Xk = f (Xk )∆t + D1/2(Xk )k i.i.d. Gaussian.

√∆t k .

Inference Problems
Given noisy observations D ≡ (y1, . . . , yN ) of hidden process X(ti)at times ti for i = 1, . . . , N .

• Estimate X(t) for 0 ≤ t ≤ T (smoothing).

• Estimate system parameters θ in drift f and diffusion D or in ratesf (X 0|X)

Can’t we just discretise ....
and run an MCMC (Gibbs) sampler ?

• Alternate between sampling paths p(x0:T |θ) and parameters p(θ|x0:T ).

• Path sampling: Paths are described by a large vectors x0:T =(x1, . . . xK ).
Sample from multivariate density
q(x0:T ) ∝ e−V (x0:T )
with
>

1 X
−1V (x0:T ) =xi+1 − xi − f (xi)∆t D (xi) xi+1 − xi − f (xi)∆t2∆t i−Xobsln p(yobs|xobs)

• Block updates using bridge proposals within Metropolis - Hastings
algorithm (Golightly, Wilkinson 2008).

20

preys1510

5510t15• Alternative: Hybrid Monte Carlo (Alexander et al 2005)

Problem with Gibbs samplers• Problem: For ∆t → 0 (complete SDE path) diffusion D and x0:T
are highly dependent:X
tXt+∆t − Xt

2→Z T
0D(Xt)dt• A Gibbs sampler will get stuck in the same (diffusion) parameters
for small ∆t.

20x(t)100-10-20

024

6810t

Ornstein Uhlenbeck - processes for σ = 1 (red), σ = 5 (green),σ = 10 (blue).

Isn’t this just some sort of hidden Markov model ?
• Yes, but we don’t know the transition probabilities.

• Can’t we still use a Forward - Backward Algorithm ?

Forward - Backward Solution• Marginal posterior at time tpt(x|{yi}N|{y }) × p({yi}ti≥t|Xt = x)i=1 ) ∝ p(x| t {zi ti<t } |}.
. {z=pF (x,t)=ψ(x,t)• Fokker - Planck equation (between observations)1∂+ ∇f − Tr(∇∇T D) pF (x, t) = 0∂t2
• Backward equation∂1+ f >∇ + Tr(D∇>∇) ψ(x, t) = 0∂t2


• & a set of jump conditions at observations.

KSP equations (Kushner ’62, Stratonovich ’60 & Pardoux ’82).

The backward equation is important !
• Posterior process is Markovian and fulfils SDE
dX(t) = g(X(t), t)dt + D1/2(X(t)) dW (t)with an effective time dependent drift g(X(t), t)g(x, t) = f (x) + D(x, t)∇ ln ψ(x, t)and ψ is the solution to the backward equation.

• Likelihood of observations P (y1, . . . , yn|X0 = x) = ψ(x, 0).

Example
Wiener process with single, noise free observation y = x(t = T ) = 0

Posterior drift g(x, t) = − T x−t for 0 < t < T .

The linear noise approximation• Similar to extended Kalman filter approach.

• Expansion in small parameter like 1/Number of molecules(Van Kampen).

• leads to linear SDE: Xt becomes approximated by a Gaussian process !

• Can be applied to backward equation → ODEs instead of PDEs !

predator population

Comparison of linear noise with MCMC

2015100510t1520

prey population

20
15
10500510t

1520

Posterior for parameters3001008020060p(α)

p(β)

40
100
20
000.020.04

α0.06

0.08

000.1

30

150

20

100p(γ)

p(δ)

1050

000.03

0.06

γ

0.09

0.12

0.15

000.01

0.005

0.005

0.015

β0.01

δ

0.015

0.02

0.02

0.025

Comments• Good approximation for small fluctuations• Does not work well for multistable systems• Gaussian approximation not optimised !

The variational approach in machine learning

Let P (X|D) =P (D|X)P (X). Then we have
P (D)P (·|D) = arg min F [Q]Qwhere the variational free energy is defined as
F [Q] = KL(Q|Pprior ) − EQ[ln P (D|·)]with the Kullback - Leibler divergence
KL(Q|P ) = EQ[lnQ] .
Pandmin F [Q] = P (D)Q

Parameter inference P (D) → p(D|θ): Approximate posterior of parameters (Lappalainen, 2000)

q(θ|D) ≈ R

e−Fθ [q] p(θ)
e−Fθ [q] p(θ) dθ.

The posterior measure
Conditional (posterior) distribution over paths X0:T given the data
D = (y1, . . . , yn)

NYPprior (X0:T )P (X0:T |D) =×p(yn|X(tn)),
P (D)
n=1

The Kullback - Leibler divergence for paths
Discretize !
Prior and posterior over paths are Markovian with transition probabilities P (x0, t0|x, tk ) and Q(x0, t0|x, tk ).

KL [QkP ] ===Z

dx0:T Q(x0:T ) ln

K−1
X Z
k=0
K−1
X Z
k=0

dx Q(x, tk )ZQ(x0:T )P (x0:T )Q(x0, tk + ∆t|x, tk )00dx Q(x , tk + ∆t|x, tk ) ln
P (x0, tk + ∆t|x, tk )dx Q(x, tk )KL [Q(·, tk + ∆t|x, tk )kP (·, tk + ∆t|x, tk )]... use short time transition probability
For diffusions

1k∆X − f (X)∆tk2
p X 0, t + ∆t|X, t ∝ exp −D2∆tas ∆t → 0,

> D −1 F .
with kF k2
=FD







and take the limit ∆t → 0
Let Q and P be measures over paths for SDEs with drifts g(X, t) andf (X, t) with same diffusion D(X). Then

1 Tdt
KL [QkP ] =2 0
Zq(x, t) is the density of X(t).

Zdx q(x, t) kg(x, t) − f (x, t)k2D

The variational problem (exact inference !)Minimise variational free energy
1 TFθ (Q) =2 0Z

Zq(x, t) kg(x, t) − f (x)k2
D dx dt −X

EQ[ln p(yi|X(ti))]i

The marginal density q fulfils the Fokker - Planck equation∂q
1= −∇g + Tr(∇∇T )D q ≡ Lg q∂t2

Variation of Lagrange function
1 Tdt dx kg(x, t) − f (x, t)k2
L =D q(x, t)
2 0 Z
X−dx q(x, ti) ln p(yi|x)]−ZZiZ TZ0dt

dx λ(x, t)∂q− Lg q∂twith respect to q and g, setting λ(x, t) = − ln ψ(x, t) yields
g(x, t) = f (x) + D(x, t)∇ ln ψ(x, t)
1∂>+ f ∇ + Tr(D∇>∇) ψ(x, t) = 0∂t2+ jump conditions at observations.

Relation to optimal control(B Kappen, E Todorov)
Consider SDE
dX(t) = (u(X(t), t) + f (X(t)))dt + D1/2(X(t)) dW (t)Adjust control u(Xt, t) so thatV =Z T0E u [Lt(Xt, ut)]

is minimal, where1Lt(x, u) = kuk2
D + U (x(t), t)
2

The Gaussian Variational Approximation

For previous applications in machine learning (see e.g. Barber & Bishop(1998), Seeger (2000), Honkela & Valpola (2005)).

The Gaussian Variational Approximation for SDEs
(Archambeau, Cornford, Opper & Shawe - Taylor, 2007)• Diffusion must be independent of X.

• Gaussian measure over paths X0:T induced by linear (approximate)posterior SDE:dX(t) = {−A(t)X + b(t)} dt + D1/2dW
Constraints are evolution eqs. for marginal mean and covariancedm= −Am + b(t)dtdS= −AS − SA> + D.
dt
→ nonlinear ODEs instead of PDEs !

Motion in double-well potential

dX = 4X(1 − X 2)dt + σ 2dW.

Variational result and comparison to MCMC
(Archambeau, Opper, Shen, Cornford, Shawe-Taylor, 2008)

Estimation of Diffusion constant σ 2q(θ|D) ≈ R

e−Fθ (Q) p(θ)
e−Fθ (Q) p(θ) dθ

1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

0.3

0.4

0.5

0.6
2!

0.7

0.8

0.9

With more observations

More dimensions
Toy weather model (Lorenz and Emanuel, 1998): x = (x1, . . . , x40) withdrifti+1i−2fi(xt) = xt− xt
xi−1− xit + θt(Vrettas, Cornford and Opper, submitted)

Large observation noise

Double well with observation noise σo = 0.6

KL divergence for Markov jump processes
Assume transition rates g(X 0|X, t) and f (X 0|X, t)

KL [Q||P ] =Z T0(dtX
xq(x, t)Xx0 :x0 6=xg(x0|x, t) lng(x0|x, t)f (x0|x, t)+ f (x0|x, t) − g(x0|x, t))Mean field approximation
Multivariate states X = (X(1), . . . , X(d))
Exact inference: Linear ODEs in S d variables
Variational approximation: Optimise in family of factorisingmeasures, i.e. of the type
Qd
Q[X0:T ] = i=1 Qi[X0:T (i)]Linear ODEs in Sd variables.
(Sanguinetti & Opper, 2008)

Simple model of autoregulatory network
2 interacting molecules: mRNA and a Protein
Number of mRNA and Protein molecules: XmRNA, XProt

X = (XmRNA, XProt)

XProt
XProt
XmRNA
XmRNA→→→→

XProt + 1
XProt − 1
XmRNA − 1XmRNA + 1::::with Ratewith Ratewith Ratewith RateγXmRNA ,δXProt ,βXmRNA ,α(1 − αcΘ(XProt − θc)) ,where Θ(x) = 1 if x ≥ 0 and Θ(x) = 0 for x < 0.

ProteinmRNA(M. Opper & G. Sanguinetti, NIPS07)
Critical parameter θc = 20 (exact) and θc = 19 (ML estimate).

Hybrid models: Inference of transcriptional regulation• Transcription factors regulate genes by binding to specific sites.

• Hard to measure transcription factor activity directly. Inference
must be based on measurement of mRNA concentration of target genes.

• Big networks: Use Clustering. Groups of genes with similar expression profiles likely regulated by same TF.

• Big networks: Latent variable models (Factor analysis).

• Small subnetworks: more detailed dynamical model (Barenco et al:)dxi
= −λixi(t) + Aif (t)dttakes sensitivity and degradation into account.
Try predictions on TF activity f (t) and learn parameters using
measurements of mRNA concentration of target genes:
yik = xi(tk ) + noise.

A priori knowledge of transcription factor activity f (t)• Discretize and assume f (t) independent for different times t (Barenco et al, 2006).

• f (t) smooth & modelled by a Gaussian process (Lawrence, Sanguinetti & Rattray, 2006)• This work: Fast switching activity µ(t) ∈ {0, 1} and µ → 1 − µ withrates γ±. (Sanguinetti, Ruttor, Archambeau, Opper 2009)• The latter model isdxi= −λixi(t) + Aiµ(t) + bidt

Telegraph Process: Xt ∈ {0, 1}x(t)1002

46t810

The role of the prior: telegraph vs GP

Variational approximation for Telegraph process• We need the variational free energy

F [Q] = KL(Q|Pprior ) − EQ[ln p(D|·)]
where Q is a Markov process over binary paths µ0:T .
• D = (y1, . . . , yn) are noisy observations of x(ti)• where x(t) is a solution todxi= −λixi(t) + Aiµ(t) + bidt
• HencenX
1p(D|µ0:T ) ∝ exp − 2(x(ti) − yi)22σ i=1• the linear ODE is solved by

x (t) = exp (−λt) x (0) +Z t0



exp (λs) [Aµ (s) + b] ds• Let f± be the prior rates and g±(t) the posterior rates. ThenZ T"#g−(t)KL QkPprior =+ f−(t) − g−(t) +dt q1(t) g−(t) lnf−(t)0"#Z Tg+(t)dt [1 − q1(t)] g+(t) ln+ f+(t) − g+(t) .
f(t)
0+hi

The total free energy

Z T"#g−(t)+ f−(t) − g−(t) +f−(t)0"#Z Tg (t)dt [1 − q1(t)] g+(t) ln ++ f+(t) − g+(t) +f+(t)0F [q] =dt q1(t) g−(t) lnn1 X
Eq2
2σ i=1

Z t0exp (−λ(t − s) [Aµ (s) + b] ds − yi2has to be minimised wrt g±(t) and q1(t) under the condition that
dq1 (t)= −(g+ + g−)q1 (t) + g+ (t)dtApproximate vs exact inferencePredictions for E. coli deprived of oxygen at t = 0based on measurements on 5 target genes at 4 time points.
(Sanguinetti et al, 2009)Multiple (2) transcription factors (Opper & Sanguinetti, 2010):

dxi= −λixi(t) + Aiµ1(t) + Aiµ1(t) + Ai12µ1(t)µ2(t) + bidt

x1

x2

0.9x3

0.7

0.8

1
0.9

0.6

0.7

0.8
0.5

0.6

0.7
0.4

0.4x0.6xx

0.5

0.5

0.3

0.3

0.4
0.2

0.2

0.3
0.1

0.1

200

300

400

500

600

700

800

900

001000100200

300

400t

500

600

700

800

900

0.1
01000

100200

300

400t

500

600

700

t

TF1

TF2

1

1

0.8

0.8

0.6

0.6

q

100

q

00

0.2

0.4

0.4

0.2

0.2

00100

200

300

400

500

t

600

700

800

900

1000

00100200

300

400

500t

600

700

800

900

1000

800

900100022Parameter inference for A21 (blue), A2 (green) and A12 (red).
!35x 1043210!1123

Prediction of activity of transcription factors FHL1 and RAP1 (Microarray data from yeast metabolic cycle)
FHL1
1
0.9
0.8
0.7

q

0.6
0.5
0.4
0.3
0.2
0.1
00100

200

300

400t

500

600

700

RAP1
10.9
0.8
0.7

q0.6
0.5
0.4
0.3
0.2
0.1
00100

200

300

400t

500

600

700

Ai1, Ai2 and Ai12 for 5 target genes.
3.5
32.5
21.5
10.5
0!0.5!1

RPL17B

RPS16B

RPL13A

RPL9A

RPL30

Comments• Method scales linearly in number of genes, number of observations& number of TFs.

• We can calculate upper and lower bounds on free energy.

• Can be generalised to more complex architectures (forward loops,autoregulatory networks).

• Can be generalised to additional white noise input (linear SDEs).

Inclusion of white noise (for single cell fluorescence measurements)

dxi = (−λixi(t) + Aiµ(t) + bi) dt + σdW
Toy data
32.5x / µ+221.5

10.500100200300

400500t

600700Opper, Ruttor and Sanguinetti, NIPS2010

8009001000

Present & Future work
• Variational path densities as proposal for MCMC• Perturbative corrections, PAC Bayesian confidence bounds• Variational approaches to multiplicative noise ?

• Realistic higher dimensional systems: Need good suboptimal variational parametrisations.