An auto-regulatory model of oscillatory gene expression

Model

Let us study the following auto-regulatory network with delay. The model is defined as follows:

\[\emptyset \xrightarrow{J_1(Y)} X,\\ Y\xrightarrow{J_2(Y)} \emptyset.\]

Note that $\emptyset \xrightarrow{J_1(Y)} X$ will trigger $X\Rightarrow Y$ after a delay time $\tau$.

According to [1], it's a simplified auto-regulatory network whereby a protein $X$ is transcribed by a gene, then it is transformed after a delay time $\tau$ into a mature protein $Y$, which binds the promoter and represses the transcription of $X$. The function $J_1(Y)$ and $J_2(Y)$ is defined as follows:

\[J_1(Y)=k_1S\frac{K^p_d}{K^p_d+Y^p},\\J_2(Y)=k_2E_T\frac{Y}{K_m+Y}.\]

In this example, we assume $k_1=k_2=S=E_T=K_d=K_m=1, p =2$ for simplicity.

Markovian part

We can define the Markovian part of the system using Catalyst

rn = @reaction_network begin
    1/(1+Y^2), 0 --> X
    1/(1+Y),   Y --> 0
end

Then we convert the reaction network to a JumpSystem

jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws=false).

We initialize the problem by setting

u0 = [0,0]
tf = 400.0
tspan = (0,tf)
τ = 20.0
dprob = DiscreteProblem(jumpsys, u0, tspan)

Non-Markovian part

We define the DelayJumpSet by

delay_trigger_affect! = function (integrator, rng)
    append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>1, 1=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

To see how to define the DelayJumpSet, we refers to this example. Thus, we can define the DelayJumpProblem by

de_chan0 = [[]]
djprob = DelayJumpProblem(jumpsys, dprob, DelayRejection(), delayjumpset, de_chan0, save_positions=(true,true)).

Solution and Visualization

Now we can solve the problem and plot two trajectories of $X$ and $Y$.

sol_1 = solve(djprob, SSAStepper(), seed = 1)
sol_2 = solve(djprob, SSAStepper(), seed = 2)

Then we simulate $10^4$ trajectories and calculate the evolution of mean value for each reactant.

using StatsBase
ens_prob = EnsembleProblem(djprob)
ens = solve(ens_prob,SSAStepper(), EnsembleThreads(), trajectories = 1e4, saveat = .1)

If we plot how the mean of $Y$ varies with respect to the mean of $X$, we will find the following oscillatory orbit in the phase diagram.

References

[1] Qingchao Jiang, Xiaoming Fu, Shifu Yan, Runlai Li, Wenli Du, Zhixing Cao, Feng Qian, Ramon Grima, "Neural network aided approximation and parameter inference of non-Markovian models of gene expression". Nature communications, (2021) 12(1), 1-12. https://doi.org/10.1038/s41467-021-22919-1