Repressilator

Repressilator model consists of a biochemical reaction network with three components in a negative feedback loop.

using Catalyst
using ModelingToolkit
using DifferentialEquations
using Plots

[ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]

repressilator = @reaction_network begin
    hillr(P₃, α, K, n), ∅ --> m₁
    hillr(P₁, α, K, n), ∅ --> m₂
    hillr(P₂, α, K, n), ∅ --> m₃
    (δ, γ), m₁ ↔ ∅
    (δ, γ), m₂ ↔ ∅
    (δ, γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    β, m₂ --> m₂ + P₂
    β, m₃ --> m₃ + P₃
    μ, P₁ --> ∅
    μ, P₂ --> ∅
    μ, P₃ --> ∅
end

\[\begin{split} \begin{align*} \require{mhchem} \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_3^{n}}} \mathrm{m_1} \\ \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_1^{n}}} \mathrm{m_2} \\ \varnothing &\xrightarrow{\frac{\alpha K^{n}}{K^{n} + P_2^{n}}} \mathrm{m_3} \\ \mathrm{m_1} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_2} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_3} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{m_1} &\xrightarrow{\beta} \mathrm{m_1} + \mathrm{P_1} \\ \mathrm{m_2} &\xrightarrow{\beta} \mathrm{m_2} + \mathrm{P_2} \\ \mathrm{m_3} &\xrightarrow{\beta} \mathrm{m_3} + \mathrm{P_3} \\ \mathrm{P_1} &\xrightarrow{\mu} \varnothing \\ \mathrm{P_2} &\xrightarrow{\mu} \varnothing \\ \mathrm{P_3} &\xrightarrow{\mu} \varnothing \end{align*} \end{split}\]

Reactions in the reaction network

reactions(repressilator)

15-element Vector{Reaction}:
 Catalyst.hillr(P₃(t), α, K, n), ∅ --> m₁
 Catalyst.hillr(P₁(t), α, K, n), ∅ --> m₂
 Catalyst.hillr(P₂(t), α, K, n), ∅ --> m₃
 δ, m₁ --> ∅
 γ, ∅ --> m₁
 δ, m₂ --> ∅
 γ, ∅ --> m₂
 δ, m₃ --> ∅
 γ, ∅ --> m₃
 β, m₁ --> m₁ + P₁
 β, m₂ --> m₂ + P₂
 β, m₃ --> m₃ + P₃
 μ, P₁ --> ∅
 μ, P₂ --> ∅
 μ, P₃ --> ∅

State variables in the reaction network

states(repressilator)

6-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 m₁(t)
 m₂(t)
 m₃(t)
 P₁(t)
 P₂(t)
 P₃(t)

Parameters in the reaction network

parameters(repressilator)

7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 α
 K
 n
 δ
 γ
 β
 μ

Convert to an ODE system

odesys = convert(ODESystem, repressilator)

\[\begin{split} \begin{align} \frac{\mathrm{d} m_1\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_3\left( t \right), \alpha, K, n \right) + \gamma - m_1\left( t \right) \delta \\ \frac{\mathrm{d} m_2\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_1\left( t \right), \alpha, K, n \right) + \gamma - m_2\left( t \right) \delta \\ \frac{\mathrm{d} m_3\left( t \right)}{\mathrm{d}t} =& \mathrm{hillr}\left( P_2\left( t \right), \alpha, K, n \right) + \gamma - m_3\left( t \right) \delta \\ \frac{\mathrm{d} P_1\left( t \right)}{\mathrm{d}t} =& - P_1\left( t \right) \mu + m_1\left( t \right) \beta \\ \frac{\mathrm{d} P_2\left( t \right)}{\mathrm{d}t} =& m_2\left( t \right) \beta - P_2\left( t \right) \mu \\ \frac{\mathrm{d} P_3\left( t \right)}{\mathrm{d}t} =& - P_3\left( t \right) \mu + m_3\left( t \right) \beta \end{align} \end{split}\]

To setup parameters (ps) and initial conditions (u₀), you can use Julia symbols to map the values.

p = [:α => 0.5, :K => 40, :n => 2, :δ => log(2) / 120, :γ => 5e-3, :β => 20 * log(2) / 120, :μ => log(2) / 60]
u₀ = [:m₁ => 0.0, :m₂ => 0.0, :m₃ => 0.0, :P₁ => 20.0, :P₂ => 0.0, :P₃ => 0.0]

6-element Vector{Pair{Symbol, Float64}}:
 :m₁ => 0.0
 :m₂ => 0.0
 :m₃ => 0.0
 :P₁ => 20.0
 :P₂ => 0.0
 :P₃ => 0.0

Or you can also use MTK symbols from the ODE system with the @unpack macro (recommended)

@unpack m₁, m₂, m₃, P₁, P₂, P₃, α, K, n, δ, γ, β, μ = repressilator
p = [α => 0.5, K => 40, n => 2, δ => log(2) / 120, γ => 5e-3, β => 20 * log(2) / 120, μ => log(2) / 60]
u₀ = [m₁ => 0.0, m₂ => 0.0, m₃ => 0.0, P₁ => 20.0, P₂ => 0.0, P₃ => 0.0]

6-element Vector{Pair{Num, Float64}}:
 m₁(t) => 0.0
 m₂(t) => 0.0
 m₃(t) => 0.0
 P₁(t) => 20.0
 P₂(t) => 0.0
 P₃(t) => 0.0

Then we can solve this reaction network

tspan = (0.0, 10000.0)
oprob = ODEProblem(odesys, u₀, tspan, p)
sol = solve(oprob)
plot(sol)

Use extracted symbols for a phase plot.

plot(sol, idxs=(P₁, P₂))

Convert to Stochastic Differential Equation (SDE) Models

Convert the very same reaction network to Chemical Langevin Equation (CLE), adding Brownian motion terms to the state variables.

tspan = (0.0, 10000.0)
sprob = SDEProblem(repressilator, u₀, tspan, p)

SDEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10000.0)
u0: 6-element Vector{Float64}:
  0.0
  0.0
  0.0
 20.0
  0.0
  0.0

solve and plot, tstops is used to specify enough points that the plot looks well-resolved.

sol = solve(sprob, LambaEulerHeun(), tstops=range(0.0, tspan[2]; length=1001))
plot(sol)

Convert to Gillespie’s stochastic simulation

Create a Gillespie stochastic simulation model from the same reaction network. The initial conditions should be integers

u₀ = [m₁ => 0, m₂ => 0, m₃ => 0, P₁ => 20, P₂ => 0, P₃ => 0]

6-element Vector{Pair{Num, Int64}}:
 m₁(t) => 0
 m₂(t) => 0
 m₃(t) => 0
 P₁(t) => 20
 P₂(t) => 0
 P₃(t) => 0

Create a discrete problem because our species are integer valued:

dprob = DiscreteProblem(repressilator, u₀, tspan, p)

DiscreteProblem with uType Vector{Int64} and tType Float64. In-place: true
timespan: (0.0, 10000.0)
u0: 6-element Vector{Int64}:
  0
  0
  0
 20
  0
  0

Create a JumpProblem, and specify Gillespie’s Direct Method as the solver.

jprob = JumpProblem(repressilator, dprob, Direct(), save_positions=(false, false))

JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 3
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 12

Solve and visualize the problem.

sol = solve(jprob, SSAStepper(), saveat=10.0)
plot(sol)