GitHub

Modelling a Periodic Event During ODE and Jump Simulations

This tutorial will describe how to simulate systems with periodic events in ODEs and jump simulations (SDEs use identical syntax). We will consider a model with a circadian rhythm, where a parameter represents the level of light. While outdoor light varies smoothly, in experimental settings a lamp is often simply turned on/off every 12 hours. Here we will model this toggling of the light using a periodic event that is triggered every 12 hours.

Modelling a circadian periodic event in an ODE simulation

We will consider a simple circadian model, consisting of a single protein ($X$), which is phosphorylated ($X \to Xᴾ$) in the presence of light ($l$). Here, the light parameter can either be $0$ (night) or $1$ (day). We can model this using a simple periodic event which switches the value of $l$ every 12 hours (here, % is the remainder operator).

using Catalyst
circadian_model = @reaction_network begin
    @discrete_events begin
        12 => [l ~ (l + 1)%2]
    end
 (kₚ*l,kᵢ), X <--> Xᴾ
end

\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[\mathtt{k_i}]{\mathtt{k_p} l} \mathrm{\mathtt{X^P}} \end{align*} \]

We can now simulate this model, observing how a 24-hour cycle is reached

using OrdinaryDiffEqDefault, Plots
u0 = [:X => 150.0, :Xᴾ => 50.0]
ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0]
tspan = (0.0, 100.0)
oprob = ODEProblem(circadian_model, u0, tspan, ps)
sol = solve(oprob)
plot(sol)
Example block output

Modelling a circadian periodic event in a jump simulation

We can define periodic events in an identical manner for jump simulations. Let's reuse our previously defined network, but now simulate it as a stochastic chemical kinetics jump process model

using JumpProcesses
u0 = [:X => 150, :Xᴾ => 50]  # define u0 as integers now
jinput = JumpInputs(circadian_model, u0, tspan, ps)
jprob = JumpProblem(jinput)
jsol = solve(jprob)
plot(jsol)
Example block output

Sometimes we might want to setup a more complicated condition than simply that our event occurs at a fixed time or with a fixed frequency. For example, suppose we want to skip the first event at time t = 12 and then have the event be periodic after that point every 12 units of time. We can do so using a more general discrete callback as follows

circadian_model = @reaction_network begin
    @discrete_events begin
        ((t % 12 == 0) & (t > 12)) => [l ~ (l + 1)%2]
    end
    (kₚ*l,kᵢ), X <--> Xᴾ
end

\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[\mathtt{k_i}]{\mathtt{k_p} l} \mathrm{\mathtt{X^P}} \end{align*} \]

Here our condition ((t % 12 == 0) & (t > 12)) determines when the event occurs, evaluating to true when t is a multiple of 12 that is also larger than 12. We now finish specifying our model

using JumpProcesses
u0 = [:X => 150, :Xᴾ => 50]
ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0]
tspan = (0.0, 100.0)
jinput = JumpInputs(circadian_model, u0, tspan, ps)
jprob = JumpProblem(jinput)

Next, if we simulate our model, we note that the events do not seem to be triggered

sol = solve(jprob)
plot(sol)
Example block output

The reason is that general discrete callbacks' conditions are only checked at the end of each simulation time step, and the probability that these will coincide with the callback's trigger times ($t = 12, 24, 36, ...$) is infinitely small. Hence, we must directly instruct our simulation to stop at these time points. We can do this using the tstops argument:

tstops = range(12.0, tspan[2]; step = 12.0)
sol = solve(jprob; tstops)
plot(sol)
Example block output

Plotting the light level

Sometimes when simulating models with periodic parameters, one would like to plot the parameter's value across the simulation. For this, there are two potential strategies. One includes creating a saving callback. The other one, which we will demonstrate here, includes turning the parameter $l$ to a variable (so that its value is recorded throughout the simulation):

circadian_model = @reaction_network begin
    @variables l(t)
    @discrete_events begin
        12 => [l ~ (l + 1)%2]
    end
 (kₚ*l,kᵢ), X <--> Xᴾ
end

Next, we simulate our model like before (but providing $l$'s value as an initial condition):

u0 = [:X => 150.0, :Xᴾ => 50.0, :l => 1.0]
ps = [:kₚ => 0.1, :kᵢ => 0.1]
oprob = ODEProblem(circadian_model, u0, (0.0, 100.0), ps)
sol = solve(oprob)

If we directly plot $l$'s value, it will be too small (compared to $X$ and $Xᴾ$ to be discernible). We instead @unpack our variables, and then plot a re-scaled version:

@unpack X, Xᴾ, l = circadian_model
plot(sol; idxs = [X, Xᴾ, 150*l], labels = ["X" "Xᴾ" "Light amplitude"])
Example block output