DiffEqBiological.jl is a domain specific language (DSL) for writing chemical reaction networks in Julia. The generated chemical reaction network model can then be translated into a variety of mathematical models which can be solved using components of the broader DifferentialEquations.jl ecosystem.
In this tutorial we'll provide an introduction to using DiffEqBiological to specify chemical reaction networks, and then to solve ODE, jump, tau-leaping and SDE models generated from them. Let's start by using the DiffEqBiological reaction_network
macro to specify a simply chemical reaction network; the well-known Repressilator.
We first import the basic packages we'll need, and use Plots.jl for making figures:
# If not already installed, first hit "]" within a Julia REPL. Then type: # add DifferentialEquations DiffEqBiological PyPlot Plots Latexify using DifferentialEquations, DiffEqBiological, Plots, Latexify pyplot(fmt=:svg);
ERROR: ArgumentError: Package PyPlot not found in current path: - Run `import Pkg; Pkg.add("PyPlot")` to install the PyPlot package.
We now construct the reaction network. The basic types of arrows and predefined rate laws one can use are discussed in detail within the DiffEqBiological Chemical Reaction Models documentation. Here we use a mix of first order, zero order and repressive Hill function rate laws. Note, $\varnothing$ corresponds to the empty state, and is used for zeroth order production and first order degradation reactions:
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 α K n δ γ β μ;
(::Main.##WeaveSandBox#253.reaction_network) (generic function with 2 metho ds)
We can use Latexify to look at the corresponding reactions and understand the generated rate laws for each reaction
latexify(repressilator; env=:chemical)
"\n\\begin{align*}\n\\require{mhchem}\n\\ce{ \\varnothing &->[\\frac{\\alph a \\cdot K^{n}}{K^{n} + P_3^{n}}] m_{1}}\\\\\n\\ce{ \\varnothing &->[\\frac {\\alpha \\cdot K^{n}}{K^{n} + P_1^{n}}] m_{2}}\\\\\n\\ce{ \\varnothing &-> [\\frac{\\alpha \\cdot K^{n}}{K^{n} + P_2^{n}}] m_{3}}\\\\\n\\ce{ m_{1} &<= >[\\delta][\\gamma] \\varnothing}\\\\\n\\ce{ m_{2} &<=>[\\delta][\\gamma] \ \varnothing}\\\\\n\\ce{ m_{3} &<=>[\\delta][\\gamma] \\varnothing}\\\\\n\\c e{ m_{1} &->[\\beta] m_{1} + P_{1}}\\\\\n\\ce{ m_{2} &->[\\beta] m_{2} + P_ {2}}\\\\\n\\ce{ m_{3} &->[\\beta] m_{3} + P_{3}}\\\\\n\\ce{ P_{1} &->[\\mu] \\varnothing}\\\\\n\\ce{ P_{2} &->[\\mu] \\varnothing}\\\\\n\\ce{ P_{3} &- >[\\mu] \\varnothing}\n\\end{align*}\n"\begin{align*} \require{mhchem} \ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + P_3^{n}}] m_{1}}\\ \ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + P_1^{n}}] m_{2}}\\ \ce{ \varnothing &->[\frac{\alpha \cdot K^{n}}{K^{n} + P_2^{n}}] m_{3}}\\ \ce{ m_{1} &<=>[\delta][\gamma] \varnothing}\\ \ce{ m_{2} &<=>[\delta][\gamma] \varnothing}\\ \ce{ m_{3} &<=>[\delta][\gamma] \varnothing}\\ \ce{ m_{1} &->[\beta] m_{1} + P_{1}}\\ \ce{ m_{2} &->[\beta] m_{2} + P_{2}}\\ \ce{ m_{3} &->[\beta] m_{3} + P_{3}}\\ \ce{ P_{1} &->[\mu] \varnothing}\\ \ce{ P_{2} &->[\mu] \varnothing}\\ \ce{ P_{3} &->[\mu] \varnothing} \end{align*}
We can also use Latexify to look at the corresponding ODE model for the chemical system
latexify(repressilator, cdot=false)
"\n\\begin{align*}\n\\frac{dm_{1}(t)}{dt} =& \\frac{\\alpha K^{n}}{K^{n} + P_3^{n}} - \\delta m_1 + \\gamma \\\\\n\\frac{dm_{2}(t)}{dt} =& \\frac{\\al pha K^{n}}{K^{n} + P_1^{n}} - \\delta m_2 + \\gamma \\\\\n\\frac{dm_{3}(t)} {dt} =& \\frac{\\alpha K^{n}}{K^{n} + P_2^{n}} - \\delta m_3 + \\gamma \\\\ \n\\frac{dP_{1}(t)}{dt} =& \\beta m_1 - \\mu P_1 \\\\\n\\frac{dP_{2}(t)}{dt } =& \\beta m_2 - \\mu P_2 \\\\\n\\frac{dP_{3}(t)}{dt} =& \\beta m_3 - \\mu P_3\n\\end{align*}\n"\begin{align*} \frac{dm_{1}(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + P_3^{n}} - \delta m_1 + \gamma \\ \frac{dm_{2}(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + P_1^{n}} - \delta m_2 + \gamma \\ \frac{dm_{3}(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + P_2^{n}} - \delta m_3 + \gamma \\ \frac{dP_{1}(t)}{dt} =& \beta m_1 - \mu P_1 \\ \frac{dP_{2}(t)}{dt} =& \beta m_2 - \mu P_2 \\ \frac{dP_{3}(t)}{dt} =& \beta m_3 - \mu P_3 \end{align*}
To solve the ODEs we need to specify the values of the parameters in the model, the initial condition, and the time interval to solve the model on. To do this it helps to know the orderings of the parameters and the species. Parameters are ordered in the same order they appear after the end
statement in the @reaction_network
macro. Species are ordered in the order they first appear within the @reaction_network
macro. We can see these orderings using the speciesmap
and paramsmap
functions:
speciesmap(repressilator)
OrderedCollections.OrderedDict{Symbol,Int64} with 6 entries: :m₁ => 1 :m₂ => 2 :m₃ => 3 :P₁ => 4 :P₂ => 5 :P₃ => 6
paramsmap(repressilator)
OrderedCollections.OrderedDict{Symbol,Int64} with 7 entries: :α => 1 :K => 2 :n => 3 :δ => 4 :γ => 5 :β => 6 :μ => 7
Knowing these orderings, we can create parameter and initial condition vectors, and setup the ODEProblem
we want to solve:
# parameters [α,K,n,δ,γ,β,μ] p = (.5, 40, 2, log(2)/120, 5e-3, 20*log(2)/120, log(2)/60) # initial condition [m₁,m₂,m₃,P₁,P₂,P₃] u₀ = [0.,0.,0.,20.,0.,0.] # time interval to solve on tspan = (0., 10000.) # create the ODEProblem we want to solve oprob = ODEProblem(repressilator, u₀, tspan, p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 10000.0) u0: [0.0, 0.0, 0.0, 20.0, 0.0, 0.0]
At this point we are all set to solve the ODEs. We can now use any ODE solver from within the DiffEq package. We'll just use the default DifferentialEquations solver for now, and then plot the solutions:
sol = solve(oprob, saveat=10.) plot(sol, fmt=:svg)
We see the well-known oscillatory behavior of the repressilator! For more on choices of ODE solvers, see the JuliaDiffEq documentation.
Let's now look at a stochastic chemical kinetics model of the repressilator, modeling it with jump processes. Here we will construct a DiffEqJump JumpProblem
that uses Gillespie's Direct
method, and then solve it to generate one realization of the jump process:
# first we redefine the initial condition to be integer valued u₀ = [0,0,0,20,0,0] # next we create a discrete problem to encode that our species are integer valued: dprob = DiscreteProblem(repressilator, u₀, tspan, p) # now we create a JumpProblem, and specify Gillespie's Direct Method as the solver: jprob = JumpProblem(dprob, Direct(), repressilator, save_positions=(false,false)) # now let's solve and plot the jump process: sol = solve(jprob, SSAStepper(), saveat=10.) plot(sol, fmt=:svg)