A toggle-switch model for gene expression

Here, we consider the toggle-switch model in Section 7.8 of Asai (2016), originated from Verd, Crombach & Jaeger (2014). See also Strasser, Theis & Marr (2012).

Toogle switches in gene expression consist of genes that mutually repress each other and exhibit two stable steady states of ON and OFF gene expression. It is a regulatory mechanism which is active during cell differentiation and is believed to act as a memory device, able to choose and maintain cell fate decisions.

The equation

We consider the following simple model as discussed in Asai (2016), of two interacting genes, $X$ and $Y$, with the concentration of their corresponding protein products at each time $t$ denoted by $X_t$ and $Y_t$. These are stochastic processes defined by the system of equations

\[ \begin{cases} \frac{\displaystyle \mathrm{d}X_t}{\displaystyle \mathrm{d} t} = \left( A_t + \frac{\displaystyle X_t^4}{\displaystyle a^4 + X_t^4}\right)\left(\frac{\displaystyle b^4}{\displaystyle b^4 + Y_t^4}\right) - \mu X_t, \\ \\ \frac{\displaystyle \mathrm{d}Y_t}{\displaystyle \mathrm{d} t} = \left( B_t + \frac{\displaystyle Y_t^4}{\displaystyle c^4 + Y_t^4}\right)\left(\frac{\displaystyle d^4}{\displaystyle d^4 + X_t^4}\right) - \nu Y_t, \\ \\ \left. X_t \right|_{t = 0} = X_0, \\ \\ \left. Y_t \right|_{t = 0} = Y_0, \end{cases}\]

where $\{A_t\}_{t\geq 0}$ and $\{B_t\}_{t\geq 0}$ are given stochastic process representing the external activation on each gene; $a$ and $c$ determine the auto-activation thresholds; $b$ and $d$ determine the thresholds for mutual repression; and $\mu$ and $\nu$ are protein decay rates. In this model, the external activation $A_t$ is a compound Poisson processes (cP), while $B_t$ is a time-dependent version of the geometric Brownian motion process (tgBm).

In the simulations below, we use the following parameters: We fix $a = c = 0.25$; $b = d = 0.4$; and $\mu = \nu = 0.75.$ The initial conditions are set to $X_0 = Y_0 = 4.0.$ The external activation $\{A_t\}_t$ is a compound Poisson process with Poisson rate $\lambda = 5.0$ and jumps uniformly distributed on $[0.0, 0.5]$. The external activation $\{B_t\}_t$ is a non-autonomous version of a geometric Brownian motion process given by

\[ \mathrm{d}B_t = (\mu_1 + \mu_2\sin(\omega t))B_t\;\mathrm{d}t + \sigma\sin(\omega t)B_t\;\mathrm{d}W_t,\]

and we choose $\mu_1 = 0.7,$ $\mu_2 = 0.3,$ $\sigma = 0.3,$ and $\omega=3\pi,$ with initial condition $A_0 = 0.2.$

We don't have an explicit solution for the equation so we just use as target for the convergence an approximate solution via Euler method at a much higher resolution.

Numerical approximation

Setting up the problem

First we load the necessary packages:

using Plots
using Measures
using Random
using LinearAlgebra
using Distributions
using RODEConvergence

Then we define the random seed:

rng = Xoshiro(123)

The evolution law:

a⁴ = c⁴ = 0.25 ^ 4
b⁴ = d⁴ = 0.4 ^ 4
μ = ν = 0.75

params = (a⁴, c⁴, b⁴, d⁴, μ,  ν)

function f!(dx, t, x, y, p)
    a⁴, c⁴, b⁴, d⁴, μ, ν = p
    α = y[1]
    β = y[2]
    x₁⁴ = x[1]^4
    x₂⁴ = x[2]^4
    dx[1] = ( α + x₁⁴ / (a⁴  + x₁⁴) ) * ( b⁴ / ( b⁴ + x₂⁴)) - μ * x[1]
    dx[2] = ( β + x₂⁴ / (c⁴  + x₂⁴) ) * ( d⁴ / ( d⁴ + x₁⁴)) - ν * x[1]
    return dx
end

The time interval:

t0, tf = 0.0, 5.0

The law for the initial condition:

x0 = 4.0
y0 = 4.0
x0law = product_distribution(Dirac(x0), Dirac(y0))

Distributions.ProductDistribution{1, 0, Tuple{Distributions.Dirac{Float64}, Distributions.Dirac{Float64}}, Distributions.Discrete, Float64}(dists=(Distributions.Dirac{Float64}(value=4.0), Distributions.Dirac{Float64}(value=4.0)), size=(2,))

The compound Poisson and the geometric Brownian motion processes, for the noisy source terms:

BM = 0.5
Bλ = 5.0
Bylaw = Uniform(0.0, BM)
Aμ1 = 0.7
Aμ2 = 0.3
Aσ = 0.3
Aω = 3π
A0 = 0.2
Aprimitive_a = t -> Aμ1 * t - Aμ2 * cos(Aω * t) / Aω
Aprimitive_b2 = t -> Aσ^2 * ( t/2 - sin(Aω * t) * cos(Aω * t) / 2Aω )
noise = ProductProcess(
    CompoundPoissonProcess(t0, tf, Bλ, Bylaw),
    HomogeneousLinearItoProcess(t0, tf, A0, Aprimitive_a, Aprimitive_b2)
)

ProductProcess{Float64, Tuple{CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}, HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}}}((CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}(0.0, 5.0, 5.0, Distributions.Uniform{Float64}(a=0.0, b=0.5)), HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}(0.0, 5.0, 0.2, Main.var"Main".var"#2#3"(), Main.var"Main".var"#5#6"())))

The resolutions for the target and approximating solutions, as well as the number of simulations for the Monte-Carlo estimate of the strong error

ntgt = 2^18
ns = 2 .^ (5:9)

5-element Vector{Int64}:
  32
  64
 128
 256
 512

nsample = ns[[1, 2, 3, 4]]

4-element Vector{Int64}:
  32
  64
 128
 256

The number of simulations for the Monte Carlo estimate is set to

m = 100

And add some information about the simulation, for the caption of the convergence figure.

info = (
    equation = "a toggle-switch model of gene regulation",
    noise = "coupled compound Poisson process and geometric Brownian motion noises",
    ic = "\$X_0 = $x0, Y_0 = $y0\$"
)

We define the target solution as the Euler approximation, which is to be computed with the target number ntgt of mesh points, and which is also the one we want to estimate the rate of convergence, in the coarser meshes defined by ns.

target = RandomEuler(length(x0law))
method = RandomEuler(length(x0law))

RandomEuler{Float64, Distributions.Multivariate}([6.94507412969136e-310, 6.94507412969215e-310])

Order of convergence

With all the parameters set up, we build the ConvergenceSuite:

suite = ConvergenceSuite(t0, tf, x0law, f!, noise, params, target, method, ntgt, ns, m)

ConvergenceSuite{Float64, Distributions.ProductDistribution{1, 0, Tuple{Distributions.Dirac{Float64}, Distributions.Dirac{Float64}}, Distributions.Discrete, Float64}, ProductProcess{Float64, Tuple{CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}, HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}}}, NTuple{6, Float64}, typeof(Main.var"Main".f!), 2, 2, RandomEuler{Float64, Distributions.Multivariate}, RandomEuler{Float64, Distributions.Multivariate}}(0.0, 5.0, Distributions.ProductDistribution{1, 0, Tuple{Distributions.Dirac{Float64}, Distributions.Dirac{Float64}}, Distributions.Discrete, Float64}(dists=(Distributions.Dirac{Float64}(value=4.0), Distributions.Dirac{Float64}(value=4.0)), size=(2,)), Main.var"Main".f!, ProductProcess{Float64, Tuple{CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}, HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}}}((CompoundPoissonProcess{Float64, Distributions.Uniform{Float64}}(0.0, 5.0, 5.0, Distributions.Uniform{Float64}(a=0.0, b=0.5)), HomogeneousLinearItoProcess{Float64, Main.var"Main".var"#2#3", Main.var"Main".var"#5#6"}(0.0, 5.0, 0.2, Main.var"Main".var"#2#3"(), Main.var"Main".var"#5#6"()))), (0.00390625, 0.00390625, 0.025600000000000005, 0.025600000000000005, 0.75, 0.75), RandomEuler{Float64, Distributions.Multivariate}([6.9451713931271e-310, 6.94521021629963e-310]), RandomEuler{Float64, Distributions.Multivariate}([6.94507412969136e-310, 6.94507412969215e-310]), 262144, [32, 64, 128, 256, 512], 100, [1, 1, 1, 1, 1], [2.55170946e-315 -9.935741340108263e-17; 2.859497335e-315 0.9160354360052033; … ; 1.448631270873832e-16 -1.436313637093215; 1.257530915583359 2.7353954553631393e-16], [2.51104501e-315 -0.15695156495366835; 3.419015454e-315 -0.15186715043709723; … ; -0.1336769169146335 0.628198574334582; -0.13378839704226084 0.631144030757221], [2.230958167e-315 5.293731102308354e180; 2.435397e-315 7.347468488156137e223; … ; 2.1723742342051675e-153 2.405388873e-315; 1.0501708915861745e224 8.74e-322])

Then we are ready to compute the errors via solve:

@time result = solve(rng, suite)

 27.967976 seconds (1.36 G allocations: 20.382 GiB, 4.03% gc time, 2.75% compilation time)

The computed strong error for each resolution in ns is stored in result.errors, and a raw LaTeX table can be displayed for inclusion in the article:

table = generate_error_table(result, suite, info)

    \begin{center}
        \begin{tabular}[htb]{|r|l|l|l|}
            \hline N & dt & error & std err \\
            \hline \hline
            32 & 0.156 & 0.645 & 0.0344 \\
            64 & 0.0781 & 0.31 & 0.0158 \\
            128 & 0.0391 & 0.15 & 0.00774 \\
            256 & 0.0195 & 0.0733 & 0.0039 \\
            512 & 0.00977 & 0.0362 & 0.00195 \\
            \hline
        \end{tabular}
    \end{center}

    \bigskip

    \caption{Mesh points (N), time steps (dt), strong error (error), and standard error (std err) of the Euler method for a toggle-switch model of gene regulation for each mesh resolution $N$, with initial condition $X_0 = 4.0, Y_0 = 4.0$ and coupled compound Poisson process and geometric Brownian motion noises, on the time interval $I = [0.0, 5.0]$, based on $M = 100$ sample paths for each fixed time step, with the target solution calculated with $262144$ points. The order of strong convergence is estimated to be $p = 1.039$, with the 95\% confidence interval $[0.965, 1.1169]$.}

The calculated order of convergence is given by result.p:

println("Order of convergence `C Δtᵖ` with p = $(round(result.p, sigdigits=2)) and 95% confidence interval ($(round(result.pmin, sigdigits=3)), $(round(result.pmax, sigdigits=3)))")

Order of convergence `C Δtᵖ` with p = 1.0 and 95% confidence interval (0.965, 1.12)

Plots

We plot the rate of convergence with the help of a plot recipe for ConvergenceResult:

plt_result = plot(result)

For the sake of illustration of the behavior of the system, we rebuild the problem with a longer time step and do a single run with it, for a single sample solution.

t0, tf = 0.0, 10.0
m = 1
ns = 2 .^ (6:9)
suite = ConvergenceSuite(t0, tf, x0law, f!, noise, params, target, method, ntgt, ns, m)
@time result = solve(rng, suite)

  0.309978 seconds (13.63 M allocations: 208.151 MiB, 7.72% gc time, 1.22% compilation time)

We visualize the pair of solutions

plt_sols = plot(suite, xshow=1, ns=nothing, label="\$X_t\$", linecolor=1)
plot!(plt_sols, suite, xshow=2, ns=nothing, label="\$Y_t\$", linecolor=2, ylabel="\$\\textrm{concentration}\$")

We also illustrate the convergence to say the expression of the X gene:

plt_suite = plot(suite, legend=:top)

We can also visualize the noises associated with this sample solution:

plt_noises = plot(suite, xshow=false, yshow=true, label=["\$A_t\$" "\$B_t\$"], linecolor=[1 2], linestyle=:auto)

We finally combine all plots into a single one, for a visual summary.

plot(plt_result, plt_sols, plt_suite, plt_noises, size=(800, 600), title=["(a) toggle-switch model" "(b) sample solution" "(c) Euler approximations" "(d) sample noises"], titlefont=10, legendfont=7)

We also combine just some of them, for the article

plt_combined = plot(plt_result, plt_suite, size=(800, 240), title=["(a) toggle-switch model" "(b) Euler approximations"], titlefont=10, layout = (1, 2), legendfont=7, bottom_margin=5mm, left_margin=5mm)

This page was generated using Literate.jl.