Product form stationary distributions of stochastic reaction - - PowerPoint PPT Presentation

Product form stationary distributions of stochastic reaction networks, and application to constrained averaging of multiscale systems

Simon Cotter

University of Manchester

19th July 2016

Simon Cotter Product form stationary dist. 0 / 35

Collaborators

Left: David Anderson (Wisconsin, US), Right: Newton Institute, University of Cambridge, UK Earlier work: Radek Erban (Oxford, UK), Kostas Zygalakis (Edinburgh, UK), Yannis Kevrekidis (Princeton, US) SLC is grateful to EPSRC for First Grant award EP/L023393/1

Simon Cotter Product form stationary dist. 1 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 1 / 35

Motivation

Simon Cotter Product form stationary dist. 2 / 35

The Random Time Change Representation

X = [X1, X2, . . . , XN] ∈ NN

0 number of each species in reactor

X(t) = X0 +

M

• m=1

Ym t αm(X(s))ds

• νm

Ym independant unit-rate Poisson processes αm(·) ≥ 0 intensities/hazard functions νm ∈ NN

0 stoichiometric vectors

Simon Cotter Product form stationary dist. 3 / 35

The Random Time Change Representation

X = [X1, X2, . . . , XN] ∈ NN

0 number of each species in reactor

X(t) = X0 +

M

• m=1

Ym t αm(X(s))ds

• νm

Ym independant unit-rate Poisson processes αm(·) ≥ 0 intensities/hazard functions νm ∈ NN

0 stoichiometric vectors

Simon Cotter Product form stationary dist. 3 / 35

The Random Time Change Representation

X = [X1, X2, . . . , XN] ∈ NN

0 number of each species in reactor

X(t) = X0 +

M

• m=1

Ym t αm(X(s))ds

• νm

Ym independant unit-rate Poisson processes αm(·) ≥ 0 intensities/hazard functions νm ∈ NN

0 stoichiometric vectors

Simon Cotter Product form stationary dist. 3 / 35

The Random Time Change Representation

X = [X1, X2, . . . , XN] ∈ NN

0 number of each species in reactor

X(t) = X0 +

M

• m=1

Ym t αm(X(s))ds

• νm

Ym independant unit-rate Poisson processes αm(·) ≥ 0 intensities/hazard functions νm ∈ NN

0 stoichiometric vectors

Simon Cotter Product form stationary dist. 3 / 35

The Random Time Change Representation

X = [X1, X2, . . . , XN] ∈ NN

0 number of each species in reactor

X(t) = X0 +

M

• m=1

Ym t αm(X(s))ds

• νm

Ym independant unit-rate Poisson processes αm(·) ≥ 0 intensities/hazard functions νm ∈ NN

0 stoichiometric vectors

Simon Cotter Product form stationary dist. 3 / 35

Gillespie SSA

[1] Calculate propensity functions αm(t), m = 1, 2, . . . , M. [2] Next reaction time is given by τ = − log (u) α0(X(t)), where α0(X(t)) =

M

• m=1

αm(X(t)). [3] Choose one j ∈ {1, . . . , M}, with probability αj/α0, and perform reaction Rj. [4] State vector X(t + τ) = X(t) + νj [5] Continue with step [1] with time t = t + τ.

Table: The Stochastic Simulation Algorithm (SSA) by Gillespie.

Simon Cotter Product form stationary dist. 4 / 35

Multiscale Systems

0.05 0.1 0.15 0.2 70 80 90 100 110 120 130 Time t Number of molecules X1 X2 X3 (X1+X2)/2

∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → X3

k5x3

− − − − → ∅ Wish to approximate by ∅

k1

− − − − → S

ˆ k4s

− − − − → X3

k5x3

− − − − → ∅

Simon Cotter Product form stationary dist. 5 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 5 / 35

The Quasi Steady-State Assumption (QSSA)

The QSSA philosophy The assumption is that the fast variables equilibrate on a timescale which is negligible with respect to the rate of change of slow species, i.e. the rates of the slow reactions can be assumed to be zero on the timescale of the fast reactions. In practice In order to approximate the effective rates of slow reactions, which are dependant on the fast variables, we only need to consider the invariant distribution of the fast reactions.

Simon Cotter Product form stationary dist. 6 / 35

The Quasi Steady-State Assumption (QSSA)

The QSSA philosophy The assumption is that the fast variables equilibrate on a timescale which is negligible with respect to the rate of change of slow species, i.e. the rates of the slow reactions can be assumed to be zero on the timescale of the fast reactions. In practice In order to approximate the effective rates of slow reactions, which are dependant on the fast variables, we only need to consider the invariant distribution of the fast reactions.

Simon Cotter Product form stationary dist. 6 / 35

QSSA: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − − − − → ← − − − − − − −

k3x2

X2

k4x2

− − − − → ∅ Effective system: ∅

k1

− − − − → S

ˆ k4s

− − − − → ∅ Fast subsystem: X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2

Simon Cotter Product form stationary dist. 7 / 35

QSSA: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − − − − → ← − − − − − − −

k3x2

X2

k4x2

− − − − → ∅ Effective system: ∅

k1

− − − − → S

ˆ k4s

− − − − → ∅ Fast subsystem: X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2

Simon Cotter Product form stationary dist. 7 / 35

QSSA: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − − − − → ← − − − − − − −

k3x2

X2

k4x2

− − − − → ∅ Effective system: ∅

k1

− − − − → S

ˆ k4s

− − − − → ∅ Fast subsystem: X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2

Simon Cotter Product form stationary dist. 7 / 35

QSSA: Simple Example

X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• k3

k2+k3 , k2 k2+k3

• steady state solution of mean field

ODE: k2λ1 = k3λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3

Simon Cotter Product form stationary dist. 8 / 35

QSSA: Simple Example

X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• k3

k2+k3 , k2 k2+k3

• steady state solution of mean field

ODE: k2λ1 = k3λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3

Simon Cotter Product form stationary dist. 8 / 35

QSSA: Simple Example

X1

k2x1

− − − − → ← − − − −

k3x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• k3

k2+k3 , k2 k2+k3

• steady state solution of mean field

ODE: k2λ1 = k3λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3

Simon Cotter Product form stationary dist. 8 / 35

QSSA: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3)s

− − − − − − − − − → ∅ Quantify error by computing error in invariant distribution of S = X1 + X2 S(∞) ∼ P k1(k2 + k3) k2k4

• (QSSA approximation)

Invariant distribution of full system given by P

• λ1 = k1

k4 , λ2 = k1(k3+k4) k2k4

• S(∞) ∼ P
• k1(k2+k3+k4)

k2k4

• (True invariant distribution)

QSSA error in Poisson intensity k1

k2

Simon Cotter Product form stationary dist. 9 / 35

QSSA: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3)s

− − − − − − − − − → ∅ Quantify error by computing error in invariant distribution of S = X1 + X2 S(∞) ∼ P k1(k2 + k3) k2k4

• (QSSA approximation)

Invariant distribution of full system given by P

• λ1 = k1

k4 , λ2 = k1(k3+k4) k2k4

• S(∞) ∼ P
• k1(k2+k3+k4)

k2k4

• (True invariant distribution)

QSSA error in Poisson intensity k1

k2

Simon Cotter Product form stationary dist. 9 / 35

QSSA: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3)s

− − − − − − − − − → ∅ Quantify error by computing error in invariant distribution of S = X1 + X2 S(∞) ∼ P k1(k2 + k3) k2k4

• (QSSA approximation)

Invariant distribution of full system given by P

• λ1 = k1

k4 , λ2 = k1(k3+k4) k2k4

• S(∞) ∼ P
• k1(k2+k3+k4)

k2k4

• (True invariant distribution)

QSSA error in Poisson intensity k1

k2

Simon Cotter Product form stationary dist. 9 / 35

QSSA: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3)s

− − − − − − − − − → ∅ Quantify error by computing error in invariant distribution of S = X1 + X2 S(∞) ∼ P k1(k2 + k3) k2k4

• (QSSA approximation)

Invariant distribution of full system given by P

• λ1 = k1

k4 , λ2 = k1(k3+k4) k2k4

• S(∞) ∼ P
• k1(k2+k3+k4)

k2k4

• (True invariant distribution)

QSSA error in Poisson intensity k1

k2

Simon Cotter Product form stationary dist. 9 / 35

QSSA: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3)s

− − − − − − − − − → ∅ Quantify error by computing error in invariant distribution of S = X1 + X2 S(∞) ∼ P k1(k2 + k3) k2k4

• (QSSA approximation)

Invariant distribution of full system given by P

• λ1 = k1

k4 , λ2 = k1(k3+k4) k2k4

• S(∞) ∼ P
• k1(k2+k3+k4)

k2k4

• (True invariant distribution)

QSSA error in Poisson intensity k1

k2

Simon Cotter Product form stationary dist. 9 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 9 / 35

The Constrained Approach

The Constrained philosophy The assumption is that the fast variables equilibrate on a timescale which is negligible with respect to the rate of change of slow species, but that the invariant distribution of the fast variables can be affected by the slow reactions. Slow stoichiometries are not in general perpendicular to the fast stoichiometries! In practice We cannot approximate P(F|S) by considering just the fast subsystem, as we need to take into account the effect of the slow reactions on the invariant distribution of the fast variables. This can be approximated by considering the dynamics of the full system, constrained to a particular value of the slow variables. NB: Assymptotically (in the limit of large timescale gaps), this approach is in agreement with the QSSA

Simon Cotter Product form stationary dist. 10 / 35

The Constrained Approach

The Constrained philosophy The assumption is that the fast variables equilibrate on a timescale which is negligible with respect to the rate of change of slow species, but that the invariant distribution of the fast variables can be affected by the slow reactions. Slow stoichiometries are not in general perpendicular to the fast stoichiometries! In practice We cannot approximate P(F|S) by considering just the fast subsystem, as we need to take into account the effect of the slow reactions on the invariant distribution of the fast variables. This can be approximated by considering the dynamics of the full system, constrained to a particular value of the slow variables. NB: Assymptotically (in the limit of large timescale gaps), this approach is in agreement with the QSSA

Simon Cotter Product form stationary dist. 10 / 35

The Constrained Approach

The Constrained philosophy The assumption is that the fast variables equilibrate on a timescale which is negligible with respect to the rate of change of slow species, but that the invariant distribution of the fast variables can be affected by the slow reactions. Slow stoichiometries are not in general perpendicular to the fast stoichiometries! In practice We cannot approximate P(F|S) by considering just the fast subsystem, as we need to take into account the effect of the slow reactions on the invariant distribution of the fast variables. This can be approximated by considering the dynamics of the full system, constrained to a particular value of the slow variables. NB: Assymptotically (in the limit of large timescale gaps), this approach is in agreement with the QSSA

Simon Cotter Product form stationary dist. 10 / 35

The Constrained Approach

R1 : ν1 → ν′

1,

R2 : ν2 → ν′

2,

. . . RM : νM → ν′

M,

Pick slow variable S =

i aijXi which is invariant with respect

to fast reactions, i.e. perpendicular to fast stoichiometries Form basis with the slow variables, along with chosen fast variables F Define constrained stoichiometric projector P : Rd → Rd P([S, F]) = [0, F]

Simon Cotter Product form stationary dist. 11 / 35

The Constrained Approach

R1 : ν1 → ν′

1,

R2 : ν2 → ν′

2,

. . . RM : νM → ν′

M,

Pick slow variable S =

i aijXi which is invariant with respect

to fast reactions, i.e. perpendicular to fast stoichiometries Form basis with the slow variables, along with chosen fast variables F Define constrained stoichiometric projector P : Rd → Rd P([S, F]) = [0, F]

Simon Cotter Product form stationary dist. 11 / 35

The Constrained Approach

R1 : ν1 → ν′

1,

R2 : ν2 → ν′

2,

. . . RM : νM → ν′

M,

Pick slow variable S =

i aijXi which is invariant with respect

to fast reactions, i.e. perpendicular to fast stoichiometries Form basis with the slow variables, along with chosen fast variables F Define constrained stoichiometric projector P : Rd → Rd P([S, F]) = [0, F]

Simon Cotter Product form stationary dist. 11 / 35

The Constrained Approach

R1 : ν1 → ν′

1,

R2 : ν2 → ν′

2,

. . . RM : νM → ν′

M,

Pick slow variable S =

i aijXi which is invariant with respect

to fast reactions, i.e. perpendicular to fast stoichiometries Form basis with the slow variables, along with chosen fast variables F Define constrained stoichiometric projector P : Rd → Rd P([S, F]) = [0, F]

Simon Cotter Product form stationary dist. 11 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

The Constrained Subsystem

R1 : ν1 → ν1 + P(ν′

1 − ν1) = ˆ

ν1, R2 : ν2 → ν2 + P(ν′

2 − ν2) = ˆ

ν2, . . . RM : νM → νM + P(ν′

M − νM) = ˆ

νM. Preserves changes to fast variables, removes changes to slow variables Multiply propensities by ✶X+P(ν′

k−νk)∈Zd 0 to preserve

non-negativity Constrains the slow variable to its starting value Invariant distribution a good approximation of P(F|S) for full system Invariant distribution of constrained system needs to be found

Simon Cotter Product form stationary dist. 12 / 35

Constrained Approach: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → ∅ First construct basis made up of slow/fast variables S = X1 + X2, F = X1

• r

F = X2 We will pick F = X2 Therefore P([X1, X2]) = [−X2, X2] Apply to all reactions to arrive at constrained subsystem

Simon Cotter Product form stationary dist. 13 / 35

Constrained Approach: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → ∅ First construct basis made up of slow/fast variables S = X1 + X2, F = X1

• r

F = X2 We will pick F = X2 Therefore P([X1, X2]) = [−X2, X2] Apply to all reactions to arrive at constrained subsystem

Simon Cotter Product form stationary dist. 13 / 35

Constrained Approach: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → ∅ First construct basis made up of slow/fast variables S = X1 + X2, F = X1

• r

F = X2 We will pick F = X2 Therefore P([X1, X2]) = [−X2, X2] Apply to all reactions to arrive at constrained subsystem

Simon Cotter Product form stationary dist. 13 / 35

Constrained Approach: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → ∅ First construct basis made up of slow/fast variables S = X1 + X2, F = X1

• r

F = X2 We will pick F = X2 Therefore P([X1, X2]) = [−X2, X2] Apply to all reactions to arrive at constrained subsystem

Simon Cotter Product form stationary dist. 13 / 35

Constrained Approach: Simple Example

Consider the system: ∅

k1

− − − − → X1

k2x1

− − − − → ← − − − −

k3x2

X2

k4x2

− − − − → ∅ First construct basis made up of slow/fast variables S = X1 + X2, F = X1

• r

F = X2 We will pick F = X2 Therefore P([X1, X2]) = [−X2, X2] Apply to all reactions to arrive at constrained subsystem

Simon Cotter Product form stationary dist. 13 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

Reaction 1: R1 : ∅

k1

− − − − → X1

P([1, 0]) = [0, 0] Constrained version of R1 has no effect on X

Reactions 2,3: R2, R3 : X1

k2x1

− − − − → ← − − − −

k3x2

X2

No change in S Projection has no effect

Reaction 4: R4 : X2

k4x2

− − − − → ∅

P([0, −1]) = [1, −1] Constrained version: R4 : X2

k4x2

− − − − − → X1

Simon Cotter Product form stationary dist. 14 / 35

Constrained Approach: Simple Example

X1

k2x1

− − − − − − − → ← − − − − − − −

(k3+k4)x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• (k3+k4)

k2+k3+k4 , k2 k2+k3+k4

• steady state solution of mean field

ODE: k2λ1 = (k3 + k4)λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3 + k4

Simon Cotter Product form stationary dist. 15 / 35

Constrained Approach: Simple Example

X1

k2x1

− − − − − − − → ← − − − − − − −

(k3+k4)x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• (k3+k4)

k2+k3+k4 , k2 k2+k3+k4

• steady state solution of mean field

ODE: k2λ1 = (k3 + k4)λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3 + k4

Simon Cotter Product form stationary dist. 15 / 35

Constrained Approach: Simple Example

X1

k2x1

− − − − − − − → ← − − − − − − −

(k3+k4)x2

X2, S = X1 + X2 Invariant distribution X2 ∼ B(S, λ2) = π(X2) [λ1, λ2] =

• (k3+k4)

k2+k3+k4 , k2 k2+k3+k4

• steady state solution of mean field

ODE: k2λ1 = (k3 + k4)λ2, λ1 + λ2 = 1 Compute expectation of the rate of reaction R4 ˆ α4 = E(α4|S) = k4E(X2|S) = k2k4S k2 + k3 + k4

Simon Cotter Product form stationary dist. 15 / 35

Constrained Approach: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3+k4)s

− − − − − − − − − − − → ∅ To compare with QSSA, compute invariant distribution S(∞) ∼ P k1(k2 + k3 + k4) k2k4

• (Constrained approximation)

Invariant distribution identical to true invariant distribution of S = X1 + X2

Simon Cotter Product form stationary dist. 16 / 35

Constrained Approach: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3+k4)s

− − − − − − − − − − − → ∅ To compare with QSSA, compute invariant distribution S(∞) ∼ P k1(k2 + k3 + k4) k2k4

• (Constrained approximation)

Invariant distribution identical to true invariant distribution of S = X1 + X2

Simon Cotter Product form stationary dist. 16 / 35

Constrained Approach: Simple Example

∅

k1

− − − − → S

k2k4/(k2+k3+k4)s

− − − − − − − − − − − → ∅ To compare with QSSA, compute invariant distribution S(∞) ∼ P k1(k2 + k3 + k4) k2k4

• (Constrained approximation)

Invariant distribution identical to true invariant distribution of S = X1 + X2

Simon Cotter Product form stationary dist. 16 / 35

Simple Example

S = X1 + X2 20 40 60 80 100 Probability Mass 0.01 0.02 0.03 0.04 0.05 0.06 0.07

Full generator QSSA generator CMA generator

SLC, “Constrained Approximation of Effective Generators for Multiscale Stochastic Reaction Networks and Application to Conditioned Path Sampling”, in review.

Simon Cotter Product form stationary dist. 17 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 17 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Complexes

S1 + S2

k

− → S3 S1 + S2 reactant complex S3 product complex C number of complexes in network 2S1

k1X1(X1−1)

− − − − → ← − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, C = 4

Simon Cotter Product form stationary dist. 18 / 35

Some Network Properties: Linkage Classes

Network can be represented uniquely by graph with C nodes 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, 2S1 S2 ∅ S1 Each connected graph is a linkage class l number of linkage classes (in this example l = 1)

Simon Cotter Product form stationary dist. 19 / 35

Some Network Properties: Linkage Classes

Network can be represented uniquely by graph with C nodes 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, 2S1 S2 ∅ S1 Each connected graph is a linkage class l number of linkage classes (in this example l = 1)

Simon Cotter Product form stationary dist. 19 / 35

Some Network Properties: Linkage Classes

Network can be represented uniquely by graph with C nodes 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, 2S1 S2 ∅ S1 Each connected graph is a linkage class l number of linkage classes (in this example l = 1)

Simon Cotter Product form stationary dist. 19 / 35

Some Network Properties: Linkage Classes

Network can be represented uniquely by graph with C nodes 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, 2S1 S2 ∅ S1 Each connected graph is a linkage class l number of linkage classes (in this example l = 1)

Simon Cotter Product form stationary dist. 19 / 35

Some Network Properties: Linkage Classes

Network can be represented uniquely by graph with C nodes 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, 2S1 S2 ∅ S1 Each connected graph is a linkage class l number of linkage classes (in this example l = 1)

Simon Cotter Product form stationary dist. 19 / 35

Some Network Properties: Stoichiometric Subspace

Denote each reaction as: νm − − − − → ν′

m

Stochiometric subspace given by: S = spanm{νm − ν′

m}

Denote s = dim(S)

Simon Cotter Product form stationary dist. 20 / 35

Some Network Properties: Stoichiometric Subspace

Denote each reaction as: νm − − − − → ν′

m

Stochiometric subspace given by: S = spanm{νm − ν′

m}

Denote s = dim(S)

Simon Cotter Product form stationary dist. 20 / 35

Some Network Properties: Stoichiometric Subspace

Denote each reaction as: νm − − − − → ν′

m

Stochiometric subspace given by: S = spanm{νm − ν′

m}

Denote s = dim(S)

Simon Cotter Product form stationary dist. 20 / 35

Some Network Properties: Deficiency Zero

Deficiency δ = C − l − s 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, δ = 4 − 1 − 2 = 1 Networks with deficiency zero have certain properties and are well studied

Simon Cotter Product form stationary dist. 21 / 35

Some Network Properties: Deficiency Zero

Deficiency δ = C − l − s 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, δ = 4 − 1 − 2 = 1 Networks with deficiency zero have certain properties and are well studied

Simon Cotter Product form stationary dist. 21 / 35

Some Network Properties: Deficiency Zero

Deficiency δ = C − l − s 2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, δ = 4 − 1 − 2 = 1 Networks with deficiency zero have certain properties and are well studied

Simon Cotter Product form stationary dist. 21 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Some Network Properties: Weak Reversibility

ODE models of zero deficiency weakly reversible networks have a unique complex balance (steady-state of the ODE) System is reversible if all reactions are reversible 1 2 3 4 5 C1 C2 C3 C4 C5

Simon Cotter Product form stationary dist. 22 / 35

Preliminaries: form of intensity functions

The theorem that follows relies on the following form of the intensity/hazard functions αm(x) = km

d

• i=1

νmi−1

• j=0

θi(xi − j), Mass action kinetics if θ is the identity

Simon Cotter Product form stationary dist. 23 / 35

Preliminaries: form of intensity functions

The theorem that follows relies on the following form of the intensity/hazard functions αm(x) = km

d

• i=1

νmi−1

• j=0

θi(xi − j), Mass action kinetics if θ is the identity

Simon Cotter Product form stationary dist. 23 / 35

Preliminaries: form of intensity functions

The theorem that follows relies on the following form of the intensity/hazard functions αm(x) = km

d

• i=1

νmi−1

• j=0

θi(xi − j), Mass action kinetics if θ is the identity

Simon Cotter Product form stationary dist. 23 / 35

Product Form Stationary Distributions

Theorem (Anderson, Craciun, Kurtz, 2010) Let {S, C, R} be a zero deficiency, weakly reversible reaction

• network. Modeled deterministically with mass action kinetics and

rate constants km the system has complex balanced equilibrium c ∈ Rd

>0. Then, the stochastically modeled system with the

specified intensity functions admits the invariant measure on each stoichiometric class π(x) ∝

d

• i=1

cxi

i

xi−1

j=0 θi(xi − j)

. The measure π can be normalized to a stationary distribution so long as it is summable.

Simon Cotter Product form stationary dist. 24 / 35

Motivating Example for New Result

2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, Multiscale system with reversible dimerisation reactions fast Consider the constrained fast subsystem with slow variable S = X1 + 2X2 2S1

k1X1(X1−1)+k3✶{X1>1}

− − − − − − − − − − − → ← − − − − − − − − − − −

(k2+k4)X2

S2, Intensity functions do not satisfy the assumptions of the existing results

Simon Cotter Product form stationary dist. 25 / 35

Motivating Example for New Result

2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, Multiscale system with reversible dimerisation reactions fast Consider the constrained fast subsystem with slow variable S = X1 + 2X2 2S1

k1X1(X1−1)+k3✶{X1>1}

− − − − − − − − − − − → ← − − − − − − − − − − −

(k2+k4)X2

S2, Intensity functions do not satisfy the assumptions of the existing results

Simon Cotter Product form stationary dist. 25 / 35

Motivating Example for New Result

2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, Multiscale system with reversible dimerisation reactions fast Consider the constrained fast subsystem with slow variable S = X1 + 2X2 2S1

k1X1(X1−1)+k3✶{X1>1}

− − − − − − − − − − − → ← − − − − − − − − − − −

(k2+k4)X2

S2, Intensity functions do not satisfy the assumptions of the existing results

Simon Cotter Product form stationary dist. 25 / 35

Motivating Example for New Result

2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, Multiscale system with reversible dimerisation reactions fast Consider the constrained fast subsystem with slow variable S = X1 + 2X2 2S1

k1X1(X1−1)+k3✶{X1>1}

− − − − − − − − − − − → ← − − − − − − − − − − −

(k2+k4)X2

S2, Intensity functions do not satisfy the assumptions of the existing results

Simon Cotter Product form stationary dist. 25 / 35

Motivating Example for New Result

2S1

k1X1(X1−1)

− − − − − − − → ← − − − − − − −

k2X2

S2, ∅

k3

− − − − → S2, S1

k4X1

− − − − → ∅, Multiscale system with reversible dimerisation reactions fast Consider the constrained fast subsystem with slow variable S = X1 + 2X2 2S1

k1X1(X1−1)+k3✶{X1>1}

− − − − − − − − − − − → ← − − − − − − − − − − −

(k2+k4)X2

S2, Intensity functions do not satisfy the assumptions of the existing results

Simon Cotter Product form stationary dist. 25 / 35

New Structural Assumption

Partition all species S = {S1, S2, . . . , SN} into S1 and S2. Si ∈ S2 if for some ηi ≥ 2 we have that νki is always some nonnegative multiple of ηi For all other Si ∈ S1, let ηi = 1. Assume that all intensity functions satisfy: αm(x) = km

d

• i=1

νmi ηi −1

• j=0

θi(Xi − jηi), where κk > 0 and θi(xi) = 0 if and only if xi ≤ αi − 1.

Simon Cotter Product form stationary dist. 26 / 35

New Structural Assumption

Partition all species S = {S1, S2, . . . , SN} into S1 and S2. Si ∈ S2 if for some ηi ≥ 2 we have that νki is always some nonnegative multiple of ηi For all other Si ∈ S1, let ηi = 1. Assume that all intensity functions satisfy: αm(x) = km

d

• i=1

νmi ηi −1

• j=0

θi(Xi − jηi), where κk > 0 and θi(xi) = 0 if and only if xi ≤ αi − 1.

Simon Cotter Product form stationary dist. 26 / 35

New Structural Assumption

Partition all species S = {S1, S2, . . . , SN} into S1 and S2. Si ∈ S2 if for some ηi ≥ 2 we have that νki is always some nonnegative multiple of ηi For all other Si ∈ S1, let ηi = 1. Assume that all intensity functions satisfy: αm(x) = km

d

• i=1

νmi ηi −1

• j=0

θi(Xi − jηi), where κk > 0 and θi(xi) = 0 if and only if xi ≤ αi − 1.

Simon Cotter Product form stationary dist. 26 / 35

New Structural Assumption

Partition all species S = {S1, S2, . . . , SN} into S1 and S2. Si ∈ S2 if for some ηi ≥ 2 we have that νki is always some nonnegative multiple of ηi For all other Si ∈ S1, let ηi = 1. Assume that all intensity functions satisfy: αm(x) = km

d

• i=1

νmi ηi −1

• j=0

θi(Xi − jηi), where κk > 0 and θi(xi) = 0 if and only if xi ≤ αi − 1.

Simon Cotter Product form stationary dist. 26 / 35

New Structural Assumption

Partition all species S = {S1, S2, . . . , SN} into S1 and S2. Si ∈ S2 if for some ηi ≥ 2 we have that νki is always some nonnegative multiple of ηi For all other Si ∈ S1, let ηi = 1. Assume that all intensity functions satisfy: αm(x) = km

d

• i=1

νmi ηi −1

• j=0

θi(Xi − jηi), where κk > 0 and θi(xi) = 0 if and only if xi ≤ αi − 1.

Simon Cotter Product form stationary dist. 26 / 35

Theorem: Product Form Stationary Distributions

Theorem (Anderson, Cotter, 2016) Let {S1 ∪ S2, C, R} be a zero deficiency, weakly reversible reaction network satisfying the above assumption. Modeled deterministically with mass action kinetics and rate constants km the system has complex balanced equilibrium c ∈ Rd

>0. Then, the

stochastically modeled system with the specified intensity functions admits the invariant measure on each stoichiometric class π(x) ∝

d

• i=1

cXi

i

⌊Xi/ηi⌋−1

j=0

θi(Xi − jηi) . The measure π can be normalized to a stationary distribution so long as it is summable.

Simon Cotter Product form stationary dist. 27 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 27 / 35

Bistable Example

R1, R2 : X2

k1x2

− − − − → ← − − − −

k2x1x2

X1 + X2, R3, R4 : ∅

k3

− − − − → ← − − − −

k4x1

X1, R5, R6 : X1 + X1

k5x1(x1−1)

− − − − → ← − − − −

k6x2

X2, R7 : X2

k7x2

− − − − → ∅. In certain parameter regimes:

Reactions R5, R6 are most frequent Invariant distribution is bimodal

Simon Cotter Product form stationary dist. 28 / 35

Bistable Example

R1, R2 : X2

k1x2

− − − − → ← − − − −

k2x1x2

X1 + X2, R3, R4 : ∅

k3

− − − − → ← − − − −

k4x1

X1, R5, R6 : X1 + X1

k5x1(x1−1)

− − − − → ← − − − −

k6x2

X2, R7 : X2

k7x2

− − − − → ∅. In certain parameter regimes:

Reactions R5, R6 are most frequent Invariant distribution is bimodal

Simon Cotter Product form stationary dist. 28 / 35

Bistable Example

R1, R2 : X2

k1x2

− − − − → ← − − − −

k2x1x2

X1 + X2, R3, R4 : ∅

k3

− − − − → ← − − − −

k4x1

X1, R5, R6 : X1 + X1

k5x1(x1−1)

− − − − → ← − − − −

k6x2

X2, R7 : X2

k7x2

− − − − → ∅. In certain parameter regimes:

Reactions R5, R6 are most frequent Invariant distribution is bimodal

Simon Cotter Product form stationary dist. 28 / 35

Bistable Example

R1, R2 : X2

k1x2

− − − − → ← − − − −

k2x1x2

X1 + X2, R3, R4 : ∅

k3

− − − − → ← − − − −

k4x1

X1, R5, R6 : X1 + X1

k5x1(x1−1)

− − − − → ← − − − −

k6x2

X2, R7 : X2

k7x2

− − − − → ∅. In certain parameter regimes:

Reactions R5, R6 are most frequent Invariant distribution is bimodal

Simon Cotter Product form stationary dist. 28 / 35

Bistable Example

Left: Log of the invariant density of the system Right: Plot of α5(X1,X2)+α6(X1,X2)

• j αj(X1,X2)

(Expectation = 69.41%)

Simon Cotter Product form stationary dist. 29 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

QSSA and Constrained Approaches: Bistable Example

Fast subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

k6x2

X2, S = X1 + 2X2 (QSSA) Slow variable S = X1 + 2X2, pick F = X2 P([X1, X2]) = [−2X2, X2] constrained stoichiometric projector Constrained subsystem: X1 + X1

k5x1(x1−1)

− − − − − − − → ← − − − − − − −

(k6+k7)x2

X2, S = X1 + 2X2 (Constrained) Use subsystems to approximate effective rates of slow reactions Stationary distributions of the subsystems can be computed via our new result

Simon Cotter Product form stationary dist. 30 / 35

Bistable Example

S = X1 + 2X2 100 200 300 400 500 Probability Mass 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018

Approximate full generator QSSA generator CMA generator

Simon Cotter Product form stationary dist. 31 / 35

Bistable Example

QSSA CMA πΩ Relative l2 error 6.347 × 10−1 1.796 × 10−2

• LH peak position

20 20 20 LH peak height 5.378 × 10−3 1.591 × 10−2 1.582 × 10−2 RH peak position 309 295 295 RH peak height 6.192 × 10−2 4.060 × 10−3 4.006 × 10−3

Simon Cotter Product form stationary dist. 32 / 35

Outline

1

Introduction

2

The Quasi Steady-State Assumption

3

The Constrained Approach

4

Product Form Stationary Distributions for Non-Mass Action Kinetics

5

Numerical Example

6

Conclusions

Simon Cotter Product form stationary dist. 32 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

Conclusions

Constrained approach better approximates P(F|S) than QSSA Corrects for the effect of the slow reactions on the distribution

• f the fast variables

Exact for monomolecular systems Can be applied iteratively, analogous to nested QSSA approaches No need for approximations or expensive simulations Leads to fast, accurate approximations of highly complex systems What constitutes the best choice of basis? Makes conditioned path sampling feasible Used in Bayesian inverse problems Constrained averaging improves reduced ODE models

Simon Cotter Product form stationary dist. 33 / 35

References

• D. Anderson & SLC, “Product-form stationary distributions for

deficiency zero networks with non-mass action kinetics” (2016) arXiv preprint arXiv:1605.07042 SLC, “Constrained approximation of effective generators for multiscale stochastic reaction networks and application to conditioned path sampling” (2015) arXiv preprint arXiv:1506.02446. SLC & R. Erban, “Error Analysis of Diffusion Approximation Methods for Multiscale Systems in Reaction Kinetics”, SIAM Journal on Scientific Computing, (2016) 38:1 144-163 10.1137/14100052X SLC, K. Zygalakis, I. Kevrekidis & R. Erban, “A constrained approach to multiscale stochastic simulation of chemically reacting systems”, The Journal of chemical physics, (2011) 135(9), 094102.

Simon Cotter Product form stationary dist. 34 / 35

Shameless Workshop Plug

Computational Challenges in Biochemical Networks: Multiscale Modelling and Inverse Problems

Tuesday 30th - Wednesday 31st August

The University of Manchester, School of Mathematjcs Alan Turing Building

Speakers: Ruth Baker, Oxford University, UK Simon Cotuer, University of Manchester, UK Radek Erban, Oxford University, UK Ramon Grima, Edinburgh University, UK Michal Komorowski, Polish Academy of Sciences, Poland Linda Petzold, UCSB, US Alexis Boukouvalas, University of Manchester, UK Chris Sherlock, Lancaster University, UK Michael Stumpf, Imperial College, UK Darren Wilkinson, Newcastle University, UK Poster session Register now at: htup://www.maths.manchester.ac.uk/news-and-events/ events/computatjonalchallenges/

Simon Cotter Product form stationary dist. 35 / 35