diff --git a/docs/src/assets/derivative_matching_cumulant.svg b/docs/src/assets/derivative_matching_cumulant.svg
index 04e74fc..5ed423d 100644
--- a/docs/src/assets/derivative_matching_cumulant.svg
+++ b/docs/src/assets/derivative_matching_cumulant.svg
@@ -1,141 +1,54 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/docs/src/tutorials/LMA_example.md b/docs/src/tutorials/LMA_example.md
index 57a0f10..75c4aa4 100644
--- a/docs/src/tutorials/LMA_example.md
+++ b/docs/src/tutorials/LMA_example.md
@@ -20,20 +20,22 @@ using Catalyst
# how the nonlinear reactions are to be transformed using LMA
rn_nonlinear = @reaction_network begin
+ @parameters σ_b σ_u ρ_b ρ_u
σ_b, g + p → 0
σ_u*(1-g), 0 ⇒ g + p
ρ_u, g → g + p
ρ_b*(1-g), 0 ⇒ p
1, p → 0
-end σ_b σ_u ρ_b ρ_u
+end
rn_linear = @reaction_network begin
+ @parameters σ_b_LMA σ_u ρ_b ρ_u
σ_b_LMA, g → 0 # typing ̄σ_b is not allowed it seems
σ_u*(1-g), 0 ⇒ g
ρ_u, g → g+p
(ρ_b*(1-g)), 0 ⇒ p
1, p → 0
-end σ_b_LMA σ_u ρ_b ρ_u
+end
```
We can now apply the LMA to find the effective parameter $\bar{σ}_b$ and generate the corresponding moment equations of the linear GRN using MomentClosure's [`linear_mapping_approximation`](@ref):
```julia
@@ -45,7 +47,7 @@ using MomentClosure
@variables g(t)
binary_vars = [speciesmap(rn_nonlinear)[g]]
-LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false)
+LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false)
display(effective_params)
```
```julia
@@ -79,7 +81,7 @@ u₀map = deterministic_IC(u₀, LMA_eqs)
oprob_LMA = ODEProblem(LMA_eqs, u₀map, tspan, p)
sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)
-plot(sol_LMA, vars=(0, [2]), label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg")
+plot(sol_LMA, idxs=[2], label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg")
```
@@ -141,11 +143,12 @@ As pointed out by the authors in [1], a more efficient way is to expand the gene
However, TaylorSeries only supports elementary function operations at the time and hence evaluating the Kummer's function $M(\cdot,\cdot,\cdot)$ requires some more work (these specialised numerics are readily available in more established scientific computing frameworks such as Mathematica but there's no fun in that). We can extend the TaylorSeries framework by constructing a function `t_pFq` that implements a recurrence relation between the Taylor coefficients for the generalized hypergeometric function `pFq` as defined in [HypergeometricFunctions.jl](https://github.com/JuliaMath/HypergeometricFunctions.jl). This can be done as follows (note that our construction is valid only for a single-variable Taylor series [`Taylor1`](https://juliadiff.org/TaylorSeries.jl/stable/api/#TaylorSeries.Taylor1)):
```julia
using TaylorSeries, HypergeometricFunctions
+using HypergeometricFunctions: pFqweniger
# please let me know if a simpler and more efficient way to do this exists!
function t_pFq(α::AbstractVector, β::AbstractVector, a::Taylor1)
order = a.order
- aux = pFq(α, β, constant_term(a))
+ aux = pFqweniger(α, β, constant_term(a))
c = Taylor1(aux, order)
iszero(order) && return c
@@ -170,13 +173,13 @@ oprob_LMA = remake(oprob_LMA, tspan=tspan)
sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)
# rebuild the symbolic expression for the effective parameter as a function of raw moments
-μ_sym = LMA_eqs.odes.states
-p_sub = Pair.(LMA_eqs.odes.ps, p)
+μ_sym = [LMA_eqs.odes.states...]
+p_sub = Dict.(Pair.(LMA_eqs.odes.ps, p)...)
avg_σ_b_sym = collect(values(effective_params))[1]
fn = build_function(substitute(avg_σ_b_sym, p_sub), μ_sym)
avg_σ_b = eval(fn)
# evaluate the time-averaged value of the effective parameter
-σ_b_avg = sum(avg_σ_b.(sol_LMA[:])) * dt / T
+@time σ_b_avg = sum(avg_σ_b.(sol_LMA.u)) * dt / T
```
We proceed with the very last steps of the LMA to obtain the probability distribution:
```julia
diff --git a/docs/src/tutorials/P53_system_example.md b/docs/src/tutorials/P53_system_example.md
index 120bd69..2d15979 100644
--- a/docs/src/tutorials/P53_system_example.md
+++ b/docs/src/tutorials/P53_system_example.md
@@ -40,7 +40,7 @@ with parameters
We begin by loading all the packages we will need
```julia
-using Catalyst, MomentClosure, OrdinaryDiffEq, DiffEqJump,
+using Catalyst, MomentClosure, OrdinaryDiffEq, JumpProcesses,
DiffEqBase.EnsembleAnalysis, Plots, Plots.PlotMeasures
```
and then build the model using Catalyst and set its parameters as follows:
@@ -48,13 +48,14 @@ and then build the model using Catalyst and set its parameters as follows:
# → for mass-actions rate
# ⇒ for non mass-actions rate
rn = @reaction_network begin
+ @parameters k₁ k₂ k₃ k₄ k₅ k₆ k₇
(k₁), 0 → x
(k₂), x → 0
(k₃*x*y/(x+k₇)), x ⇒ 0
(k₄*x), 0 ⇒ y₀
(k₅), y₀ → y
(k₆), y → 0
-end k₁ k₂ k₃ k₄ k₅ k₆ k₇
+end
# parameters [k₁, k₂, k₃, k₄, k₅, k₆, k₇]
p = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01]
@@ -106,7 +107,7 @@ h2 = histogram(data[2], normalize=true, xlabel="y₀", ylabel="P(y₀)")
h3 = histogram(data[3], normalize=true, xlabel="y", ylabel="P(y)")
using Plots.PlotMeasures
plot(h1, h2, h3, legend=false, layout=(1,3), size = (1050, 250),
- left_margin = 5mm, bottom_margin = 7mm, guidefontsize=10)
+ left_margin = 5PlotMeasures.mm, bottom_margin = 7PlotMeasures.mm, guidefontsize=10)
```
@@ -120,7 +121,7 @@ closures = ["normal", "log-normal", "gamma"]
plts = [plot() for i in 1:length(closures)]
for q in 3:6
- eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false)
+ eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)
for (closure, plt) in zip(closures, plts)
closed_eqs = moment_closure(eqs, closure)
@@ -128,7 +129,7 @@ for q in 3:6
oprob = ODEProblem(closed_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
- plt = plot!(plt, sol, vars=(0, 1), lw=3, label = "q = "*string(q))
+ plt = plot!(plt, sol, idxs=[1], lw=3, label = "q = "*string(q))
end
end
@@ -139,13 +140,13 @@ end
```
Normal closure:
```julia
-plot(plts[1], size=(750, 450), leftmargin=2mm)
+plot(plts[1], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
Log-normal closure:
```julia
-plot(plts[2], size=(750, 450), leftmargin=2mm)
+plot(plts[2], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
@@ -157,7 +158,7 @@ plot(plts[2], xlims=(0., 50.), lw=3)
Gamma closure:
```julia
-plot(plts[3], size=(750, 450), leftmargin=2mm)
+plot(plts[3], size=(750, 450), leftmargin=2PlotMeasures.mm)
```
@@ -168,7 +169,7 @@ We can also plot the variance predictions:
# rerunning the same calculations as they are reasonably fast
plt = plot()
for q in [4,6]
- eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false)
+ eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)
for closure in closures
closed_eqs = moment_closure(eqs, closure)
@@ -177,7 +178,7 @@ for q in [4,6]
sol = solve(oprob, Tsit5(), saveat=0.1)
# index of M₂₀₀ can be checked with `u₀map` or `closed_eqs.odes.states`
- plt = plot!(plt, sol, vars=(0, 4), lw=3, label = closure*" q = "*string(q))
+ plt = plot!(plt, sol, idxs=[4], lw=3, label = closure*" q = "*string(q))
end
end
@@ -200,7 +201,7 @@ q_vals = [4, 6]
for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars)
- eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaw=false)
+ eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaws=false)
for closure in closures
closed_eqs = moment_closure(eqs, closure)
@@ -209,8 +210,8 @@ for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars)
oprob = ODEProblem(closed_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
- plt_m = plot!(plt_m, sol, vars=(0, 1), label = closure)
- plt_v = plot!(plt_v, sol, vars=(0, 4), label = closure)
+ plt_m = plot!(plt_m, sol, idxs=[1], label = closure)
+ plt_v = plot!(plt_v, sol, idxs=[4], label = closure)
end
@@ -270,7 +271,7 @@ Finally, we can check whether better estimates can be obtained using even higher
plt = plot()
closures = ["zero", "normal", "log-normal", "gamma"]
-eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaw=false)
+eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaws=false)
# faster to store than recompute in case we want to try different solvers/params
oprobs = Dict()
@@ -281,7 +282,7 @@ for closure in closures
oprobs[closure] = ODEProblem(closed_eqs, u₀map, tspan, p)
sol = solve(oprobs[closure], Tsit5(), saveat=0.1)
- plt = plot!(plt, sol, vars=(0, 1), label = closure)
+ plt = plot!(plt, sol, idxs=[1], label = closure)
end
plt = plot!(plt, xlabel = "Time [h]", ylabel = "Mean p53 molecule number")
diff --git a/docs/src/tutorials/common_issues.md b/docs/src/tutorials/common_issues.md
index ef7452e..1811970 100644
--- a/docs/src/tutorials/common_issues.md
+++ b/docs/src/tutorials/common_issues.md
@@ -7,16 +7,17 @@ We first redefine the system and its parameters for completeness:
using MomentClosure, Catalyst, OrdinaryDiffEq, Plots
rn = @reaction_network begin
+ @parameters c₁ c₂ c₃ c₄ Ω
(c₁/Ω^2), 2X + Y → 3X
(c₂), X → Y
(c₃*Ω, c₄), 0 ↔ X
-end c₁ c₂ c₃ c₄ Ω
+end
p = [0.9, 2, 1, 1, 100]
u₀ = [1, 1]
tspan = (0., 100.)
-raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false)
+raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)
```
As we have seen earlier, second-order moment expansion using normal closure approximates the true system dynamics sufficiently accurately but it's interesting to see how other closures compare. Let's try applying zero closure:
```julia
@@ -26,7 +27,7 @@ u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
@@ -40,7 +41,7 @@ u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2, legend=:bottomright)
+plot(sol, idxs=[1,2], lw=2, legend=:bottomright)
```
@@ -48,66 +49,66 @@ We observe sustained oscillatory behaviour instead of the expected damped oscill
Normal closure is also quite fragile. This can be seen by simply including the combinatorial scaling of the mass-action propensity functions with `combinatoric_ratelaw=true` which leads to unphysical sustained oscillatory trajectories:
```julia
-raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=true)
+raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true)
closed_raw_eqs = moment_closure(raw_eqs, "normal")
u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
Nevertheless, this can be improved upon by increasing the order of moment expansion:
```julia
-raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaw=true)
+raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true)
closed_raw_eqs = moment_closure(raw_eqs, "normal")
u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2, legend=:bottomright)
+plot(sol, idxs=[1,2], lw=2, legend=:bottomright)
```
Some dampening in the system is now visible. Increasing the expansion order to `4` finally leads to physically sensible results:
```julia
-raw_eqs = generate_raw_moment_eqs(rn, 4, combinatoric_ratelaw=true)
+raw_eqs = generate_raw_moment_eqs(rn, 4, combinatoric_ratelaws=true)
closed_raw_eqs = moment_closure(raw_eqs, "normal")
u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Tsit5(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
For dessert, we consider unphysical divergent trajectories—a frequent problem with moment equations [3]. A good example is the second-order moment expansion including the combinatorial scaling of propensities with log-normal closure applied:
```julia
-raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=true)
+raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true)
closed_raw_eqs = moment_closure(raw_eqs, "log-normal")
u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Rodas4P(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
In contrast to normal closure, increasing the expansion order makes the problem worse:
```julia
-raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaw=true)
+raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true)
closed_raw_eqs = moment_closure(raw_eqs, "log-normal")
u₀map = deterministic_IC(u₀, closed_raw_eqs)
oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
sol = solve(oprob, Rodas4P(), saveat=0.1)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
```julia
┌ Warning: Interrupted. Larger maxiters is needed.
diff --git a/docs/src/tutorials/derivative_matching_example.md b/docs/src/tutorials/derivative_matching_example.md
index cb5704c..57eace8 100644
--- a/docs/src/tutorials/derivative_matching_example.md
+++ b/docs/src/tutorials/derivative_matching_example.md
@@ -12,9 +12,10 @@ The reaction network and its parameters can be defined as follows:
using Catalyst
rn = @reaction_network begin
+ @parameters c₁ c₂
(c₁), x₁ → 2x₁+x₂
(c₂), x₁+x₂ → x₂
-end c₁ c₂
+end
# parameter values
p = [1.0, 1.0]
@@ -133,7 +134,7 @@ true
```
The last ingredient we need for a proper comparison between the second and third order moment expansions is a reference value predicted by the SSA. We can simulate $10^5$ SSA trajectories as follows:
```julia
-using DiffEqJump
+using JumpProcesses
dprob = DiscreteProblem(rn, u0, tspan, p)
jprob = JumpProblem(rn, dprob, Direct(), save_positions=(false, false))
diff --git a/docs/src/tutorials/geometric_reactions+conditional_closures.md b/docs/src/tutorials/geometric_reactions+conditional_closures.md
index 31f2134..3e26e7a 100644
--- a/docs/src/tutorials/geometric_reactions+conditional_closures.md
+++ b/docs/src/tutorials/geometric_reactions+conditional_closures.md
@@ -28,11 +28,12 @@ using MomentClosure, Catalyst, Distributions, JumpProcesses, DiffEqBase, Ordinar
m = rand(Distributions.Geometric(1/(1+b)))
rn = @reaction_network begin
+ @parameters k_on k_off k_p γ_p
k_on*(1-g), 0 --> g # G* -> G
k_off*P^2, g --> 0 # G -> G*
k_p, g --> g + $m*P # G -> G + mP, m ~ Geometric(p)
γ_p, P --> 0 # P -> ∅
-end k_on k_off k_p γ_p
+end
```
We can now generate the raw moment equations up to third order:
```julia
@@ -40,62 +41,18 @@ eqs = generate_raw_moment_eqs(rn, 3)
latexify(eqs)
```
```math
-
\begin{align*}
-\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\
+\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{off} \mu_{1 2} - k_{on} \mu_{1 0} \\
\frac{d\mu_{0 1}}{dt} =& b k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\
-\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{on} \mu_{1 0} + k_{off} \mu_{1 2} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\
-\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{2 0} - k_{off} \mu_{1 3} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} \\
-\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} + 2 b k_{p} \mu_{1 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} - k_{p} \mu_{1 0}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + 2 b k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} \\
-\frac{d\mu_{3 0}}{dt} =& k_{on} + 2 k_{on} \mu_{1 0} + 3 k_{off} \mu_{2 2} - k_{off} \mu_{1 2} - 3 k_{on} \mu_{3 0} - 3 k_{off} \mu_{3 2} \\
-\frac{d\mu_{2 1}}{dt} =& k_{off} \mu_{1 3} + k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + b k_{p} \mu_{3 0} - 2 k_{off} \mu_{2 3} - \gamma_{p} \mu_{2 1} - 2 k_{on} \mu_{2 1} \\
-\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{2 0} + 2 b k_{p} \mu_{2 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{2 0} - k_{p} \mu_{2 0}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + 2 b k_{p} \mu_{2 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\
-\frac{d\mu_{0 3}}{dt} =& \frac{3 k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 1} + 6 b k_{p} \mu_{1 1} + 3 b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 1} - 3 k_{p} \mu_{1 1}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{3} \mu_{1 0} + 6 b k_{p} \mu_{1 0} + 7 k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{3} \mu_{1 0} + 7 b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} - 6 k_{p} \mu_{1 0}}{b \left( \frac{1}{1 + b} \right)^{3} + \left( \frac{1}{1 + b} \right)^{3}} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 2} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3}
-\end{align*}
-```
-Note that certain symbolic expressions have not been fully simplified: this happens as [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) used internally is not yet able to simplify more complicated fractions to the extent that is required here. For reference, the simplified moment equations take the form:
-```math
-\begin{align*}
-\frac{d\mu{_{10}}}{dt} =& k_{on} - k_{off} \mu{_{12}} - k_{on} \mu{_{10}} \\
-\frac{d\mu{_{01}}}{dt} =& b k_{p} \mu{_{10}} - \gamma_{p} \mu{_{01}} \\
-\frac{d\mu{_{20}}}{dt} =& k_{on} + k_{off} \mu{_{12}} + k_{on} \mu{_{10}} - 2 k_{off} \mu{_{22}} - 2 k_{on} \mu{_{20}} \\
-\frac{d\mu{_{11}}}{dt} =& k_{on} \mu{_{01}} + b k_{p} \mu{_{20}} - k_{off} \mu{_{13}} - k_{on} \mu{_{11}} - \gamma_{p} \mu{_{11}} \\
-\frac{d\mu{_{02}}}{dt} =& \gamma_{p} \mu{_{01}} + b k_{p} \mu{_{10}} + 2 b k_{p} \mu{_{11}} + 2 k_{p} \mu{_{10}} b^{2} - 2 \gamma_{p} \mu{_{02}} \\
-\frac{d\mu{_{30}}}{dt} =& k_{on} + 3 k_{off} \mu{_{22}} + 2 k_{on} \mu{_{10}} - k_{off} \mu{_{12}} - 3 k_{off} \mu{_{32}} - 3 k_{on} \mu{_{30}} \\
-\frac{d\mu{_{21}}}{dt} =& k_{off} \mu{_{13}} + k_{on} \mu{_{01}} + k_{on} \mu{_{11}} + b k_{p} \mu{_{30}} - 2 k_{off} \mu{_{23}} - 2 k_{on} \mu{_{21}} - \gamma_{p} \mu{_{21}} \\
-\frac{d\mu{_{12}}}{dt} =& k_{on} \mu{_{02}} + \gamma_{p} \mu{_{11}} + b k_{p} \mu{_{20}} + 2 b k_{p} \mu{_{21}} + 2 k_{p} \mu{_{20}} b^{2} - k_{off} \mu{_{14}} - k_{on} \mu{_{12}} - 2 \gamma_{p} \mu{_{12}} \\
-\frac{d\mu{_{03}}}{dt} =& b k_{p} \mu{_{10}} + 3 \gamma_{p} \mu{_{02}} + 3 b k_{p} \mu{_{11}} + 3 b k_{p} \mu{_{12}} + 6 k_{p} \mu{_{10}} b^{2} + 6 k_{p} \mu{_{10}} b^{3} + 6 k_{p} \mu{_{11}} b^{2} - \gamma_{p} \mu{_{01}} - 3 \gamma_{p} \mu{_{03}}
-\end{align*}
-```
-The fractions appeared in our equations as we defined the success probability $p$ of the geometric distribution as $p = \frac{1}{1+b}$. This could have been avoided by leaving $p$ in the symbolic expressions and substituting its numerical value expressed in terms of the mean $b$ later on. We proceed to do that and redefine the reaction system accordingly, making the moment equations easier to handle:
-```julia
-@parameters p
-m = rand(Distributions.Geometric(p))
-
-rn = @reaction_network begin
- k_on*(1-g), 0 --> g
- k_off*P^2, g --> 0*
- k_p, g --> g + $m*P
- γ_p, P --> 0
-end k_on k_off k_p γ_p
-
-eqs = generate_raw_moment_eqs(rn, 3)
-latexify(eqs)
-```
-```math
-\begin{align*}
-\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\
-\frac{d\mu_{0 1}}{dt} =& k_{p} p^{-1} \mu_{1 0} - k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\
-\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{on} \mu_{1 0} + k_{off} \mu_{1 2} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\
-\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + k_{p} p^{-1} \mu_{2 0} - k_{off} \mu_{1 3} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{p} \mu_{2 0} \\
-\frac{d\mu_{0 2}}{dt} =& k_{p} \mu_{1 0} + \gamma_{p} \mu_{0 1} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} - 3 k_{p} p^{-1} \mu_{1 0} \\
-\frac{d\mu_{3 0}}{dt} =& k_{on} + 2 k_{on} \mu_{1 0} + 3 k_{off} \mu_{2 2} - k_{off} \mu_{1 2} - 3 k_{on} \mu_{3 0} - 3 k_{off} \mu_{3 2} \\
-\frac{d\mu_{2 1}}{dt} =& k_{off} \mu_{1 3} + k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + k_{p} p^{-1} \mu_{3 0} - 2 k_{off} \mu_{2 3} - \gamma_{p} \mu_{2 1} - k_{p} \mu_{3 0} - 2 k_{on} \mu_{2 1} \\
-\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + k_{p} \mu_{2 0} + 2 k_{p} p^{-1} \mu_{2 1} + 2 k_{p} p^{-2} \mu_{2 0} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} - 2 k_{p} \mu_{2 1} - 3 k_{p} p^{-1} \mu_{2 0} \\
-\frac{d\mu_{0 3}}{dt} =& 3 k_{p} \mu_{1 1} + 3 \gamma_{p} \mu_{0 2} + 7 k_{p} p^{-1} \mu_{1 0} + 3 k_{p} p^{-1} \mu_{1 2} + 6 k_{p} p^{-3} \mu_{1 0} + 6 k_{p} p^{-2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} - 3 k_{p} \mu_{1 2} - 3 \gamma_{p} \mu_{0 3} - 9 k_{p} p^{-1} \mu_{1 1} - 12 k_{p} p^{-2} \mu_{1 0}
+\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{off} \mu_{1 2} + k_{on} \mu_{1 0} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\
+\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{2 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{off} \mu_{1 3} \\
+\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - 2 \gamma_{p} \mu_{0 2} \\
+\frac{d\mu_{3 0}}{dt} =& k_{on} + 3 k_{off} \mu_{2 2} + 2 k_{on} \mu_{1 0} - k_{off} \mu_{1 2} - 3 k_{off} \mu_{3 2} - 3 k_{on} \mu_{3 0} \\
+\frac{d\mu_{2 1}}{dt} =& k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + k_{off} \mu_{1 3} + b k_{p} \mu_{3 0} - 2 k_{on} \mu_{2 1} - \gamma_{p} \mu_{2 1} - 2 k_{off} \mu_{2 3} \\
+\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + b k_{p} \mu_{2 0} + 2 k_{p} b^{2} \mu_{2 0} + 2 b k_{p} \mu_{2 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\
+\frac{d\mu_{0 3}}{dt} =& b k_{p} \mu_{1 0} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 1} + 3 b k_{p} \mu_{1 2} + 6 k_{p} b^{2} \mu_{1 0} + 6 k_{p} b^{3} \mu_{1 0} + 6 k_{p} b^{2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3}
\end{align*}
```
-
A lot of information in this system of ODEs is redundant as the gene state is a Bernoulli variable that (in our case) has the following properties:
```math
\begin{align*}
@@ -113,12 +70,12 @@ latexify(clean_eqs)
```
```math
\begin{align*}
-\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\
-\frac{d\mu_{0 1}}{dt} =& k_{p} p^{-1} \mu_{1 0} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} \\
-\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + k_{p} p^{-1} \mu_{1 0} - k_{off} \mu_{1 3} - k_{p} \mu_{1 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} \\
-\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + k_{p} \mu_{1 0} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} - 3 k_{p} p^{-1} \mu_{1 0} \\
-\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + k_{p} \mu_{1 0} + \gamma_{p} \mu_{1 1} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} - 3 k_{p} p^{-1} \mu_{1 0} \\
-\frac{d\mu_{0 3}}{dt} =& 3 k_{p} \mu_{1 1} + 3 \gamma_{p} \mu_{0 2} + 7 k_{p} p^{-1} \mu_{1 0} + 3 k_{p} p^{-1} \mu_{1 2} + 6 k_{p} p^{-3} \mu_{1 0} + 6 k_{p} p^{-2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} - 3 k_{p} \mu_{1 2} - 3 \gamma_{p} \mu_{0 3} - 9 k_{p} p^{-1} \mu_{1 1} - 12 k_{p} p^{-2} \mu_{1 0}
+\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{off} \mu_{1 2} - k_{on} \mu_{1 0} \\
+\frac{d\mu_{0 1}}{dt} =& b k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\
+\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{1 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{off} \mu_{1 3} \\
+\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - 2 \gamma_{p} \mu_{0 2} \\
+\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\
+\frac{d\mu_{0 3}}{dt} =& b k_{p} \mu_{1 0} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 1} + 3 b k_{p} \mu_{1 2} + 6 k_{p} b^{2} \mu_{1 0} + 6 k_{p} b^{3} \mu_{1 0} + 6 k_{p} b^{2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3}
\end{align*}
```
The system of ODEs is now much simpler and we can see that there are two higher-order moments we need to truncate: $\mu_{13}$ and $μ_{14}$. We consider normal, derivative matching, conditional gaussian and conditional derivative matching closures to see how well they compare. First we apply different closures and print out the corresponding higher-order moment expressions in order to check that our results are consistent with those published in [1].
@@ -143,8 +100,8 @@ latexify(dm_eqs, :closure)
```
```math
\begin{align*}
-\mu{_{13}} =& \mu{_{03}} \mu{_{10}} \mu{_{01}}^{3} \mu{_{02}}^{-3} \mu{_{11}}^{-3} \mu{_{12}}^{3} \\
-\mu{_{14}} =& \mu{_{04}} \mu{_{01}}^{-4} \mu{_{02}}^{6} \mu{_{03}}^{-4} \mu{_{10}}^{-1} \mu{_{11}}^{4} \mu{_{12}}^{-6} \mu{_{13}}^{4}
+\mu_{1 3} =& \frac{\mu_{0 1}^{3} \mu_{1 2}^{3} \mu_{0 3} \mu_{1 0}}{\mu_{0 2}^{3} \mu_{1 1}^{3}} \\
+\mu_{1 4} =& \frac{\mu_{0 2}^{6} \mu_{1 1}^{4} \mu_{1 3}^{4} \mu_{0 4}}{\mu_{0 1}^{4} \mu_{0 3}^{4} \mu_{1 2}^{6} \mu_{1 0}}
\end{align*}
```
For conditional gaussian closure:
@@ -154,8 +111,8 @@ latexify(cond_gaussian_eqs, :closure)
```
```math
\begin{align*}
-\mu{_{13}} =& 3 \mu{_{11}} \mu{_{12}} \mu{_{10}}^{-1} - 2 \mu{_{10}}^{-2} \mu{_{11}}^{3} \\
-\mu{_{14}} =& 3 \mu{_{10}}^{-1} \mu{_{12}}^{2} + 6 \mu{_{10}}^{-3} \mu{_{11}}^{4} + 4 \mu{_{11}} \mu{_{13}} \mu{_{10}}^{-1} - 12 \mu{_{12}} \mu{_{10}}^{-2} \mu{_{11}}^{2}
+\mu_{1 3} =& \frac{-2 \mu_{1 1}^{3}}{\mu_{1 0}^{2}} + \frac{3 \mu_{1 1} \mu_{1 2}}{\mu_{1 0}} \\
+\mu_{1 4} =& \frac{3 \mu_{1 2}^{2}}{\mu_{1 0}} + \frac{6 \mu_{1 1}^{4}}{\mu_{1 0}^{3}} + \frac{-12 \mu_{1 1}^{2} \mu_{1 2}}{\mu_{1 0}^{2}} + \frac{4 \mu_{1 1} \mu_{1 3}}{\mu_{1 0}}
\end{align*}
```
And, finally, for conditional derivative matching:
@@ -165,15 +122,14 @@ latexify(cond_dm_eqs, :closure)
```
```math
\begin{align*}
-\mu{_{13}} =& \mu{_{10}} \mu{_{11}}^{-3} \mu{_{12}}^{3} \\
-\mu{_{14}} =& \mu{_{10}}^{-1} \mu{_{11}}^{4} \mu{_{12}}^{-6} \mu{_{13}}^{4}
+\mu_{1 3} =& \frac{\mu_{0 1}^{3} \mu_{1 2}^{3} \mu_{0 3} \mu_{1 0}}{\mu_{0 2}^{3} \mu_{1 1}^{3}} \\
+\mu_{1 4} =& \frac{\mu_{0 2}^{6} \mu_{1 1}^{4} \mu_{1 3}^{4} \mu_{0 4}}{\mu_{0 1}^{4} \mu_{0 3}^{4} \mu_{1 2}^{6} \mu_{1 0}}
\end{align*}
```
All these results are consistent with [1], reassuring that the model and closures are implemented as intended. Finally, we can proceed to solve the resulting ODEs and compare the resulting means and standard deviations to the SSA. Following Soltani et al. [1], we define all model parameters as:
```julia
mean_p = 200
mean_b = 70
-p_val = 1/(1+mean_b) # note we use p = 1/(1+b)
γ_p_val = 1
k_off_val = 0.001
k_on_val = 0.05
@@ -184,7 +140,7 @@ symmap = [:k_on => k_on_val,
:k_off => k_off_val,
:k_p => k_p_val,
:γ_p => γ_p_val,
- :p => p_val]
+ :b => mean_b]
pmap = symmap_to_varmap(rn, symmap)
# initial gene state and protein number, order [g, p]
@@ -234,7 +190,7 @@ for closure in ["normal", "derivative matching",
# μ₀₁ is 2nd and μ₀₂ is 4th element in sol
# can check the order with `closed_eqs.odes.states`
- plt_m = plot!(plt_m, sol, vars=(0, 2), label=closure)
+ plt_m = plot!(plt_m, sol, idxs=[2], label=closure)
plt_std = plot!(plt_std, sol.t, sqrt.(sol[4, :] .- sol[2, :].^2), label=closure)
end
@@ -267,6 +223,7 @@ l = rand(Distributions.Geometric(1/(1+b_y)))
# g_y - gene state of Y protein producing gene
# x, y - proteins X and Y
rn = @reaction_network begin
+ @parameters kx_on kx_off ky_on ky_off k_x γ_x k_y γ_y
kx_on*(1-g_x)*y, 0 --> g_x # 0 -> g_x
kx_off, g_x --> 0 # g_x -> 0
ky_on*(1-g_y), 0 --> g_y # 0 -> g_y
@@ -275,7 +232,7 @@ rn = @reaction_network begin
γ_x, x --> 0 # x -> 0
k_y*g_y, 0 --> $l*y # 0 -> ly, l ~ Geometric(mean_b_y)
γ_y, y --> 0 # y -> 0
-end kx_on kx_off ky_on ky_off k_x γ_x k_y γ_y
+end
# both g_x and g_y are Bernoulli random variables
binary_vars = [1, 2]
@@ -351,7 +308,7 @@ for closure in ["derivative matching", "conditional derivative matching"]
sol = solve(oprob, Tsit5(), saveat=0.1)
# μ₀₀₀₁ is the 4th and μ₀₀₀₂ is the 12th element in sol (can check with closed_eqs.odes.states)
- plt_m = plot!(plt_m, sol, vars=(0, 4), label=closure)
+ plt_m = plot!(plt_m, sol, idxs=[4], label=closure)
plt_std = plot!(plt_std, sol.t, sqrt.(sol[12, :] .- sol[4, :].^2), label=closure)
end
diff --git a/docs/src/tutorials/parameter_estimation_SDE.md b/docs/src/tutorials/parameter_estimation_SDE.md
index 1c00b3e..1988cbe 100644
--- a/docs/src/tutorials/parameter_estimation_SDE.md
+++ b/docs/src/tutorials/parameter_estimation_SDE.md
@@ -111,7 +111,7 @@ closed_moment_prob = ODEProblem(LV_moments, u0map, (0.0, Tf), zeros(5))
function obj_MCA(p)
prob = remake(closed_moment_prob; p=p)
sol = solve(prob, Tsit5(), saveat = t_data)
- if sol.retcode == :Success
+ if sol.retcode == ReturnCode.Success
obj = sum(norm(sol.u[i][1:2] - means[i])^2 for i in 1:length(t_data))
obj += 1e4*sum((sol.u[i][3] - sol.u[i][1]^2 - vars[i][1])^2 for i in 1:length(t_data))
obj += 1e4*sum((sol.u[i][5] - sol.u[i][2]^2 - vars[i][2])^2 for i in 1:length(t_data))
diff --git a/docs/src/tutorials/time-dependent_propensities.md b/docs/src/tutorials/time-dependent_propensities.md
index 58f8ee0..5e360e5 100644
--- a/docs/src/tutorials/time-dependent_propensities.md
+++ b/docs/src/tutorials/time-dependent_propensities.md
@@ -20,12 +20,13 @@ where $c_2^0$ is a fixed value and the sinusoidal modulation with frequency $\om
using Catalyst, MomentClosure, Latexify
rn = @reaction_network begin
+ @parameters c₁ c₂ c₃ c₄ Ω ω τ
(c₁/Ω^2), 2X + Y → 3X
(c₂*(1+0.5*sin(ω*(t<τ)*t))), X → Y
(c₃*Ω, c₄), 0 ↔ X
-end c₁ c₂ c₃ c₄ Ω ω τ
+end
-raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false)
+raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)
latexify(raw_eqs)
```
```math
@@ -60,13 +61,13 @@ oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)
# solve using Tsit5 solver
sol = solve(oprob, Tsit5(), saveat=0.2)
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
It would be great to compare our results to the true dynamics. Using DifferentialEquations, we can run a modified SSA taking into account the time-dependent propensity functions ([VariableRateJumps](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/#Adding-a-VariableRateJump)). This requires some care and is done by combining an `ODEProblem` with a `JumpProblem`:
```julia
-using DiffEqJump
+using JumpProcesses
# convert ReactionSystem into JumpSystem
jsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)
diff --git a/docs/src/tutorials/using_momentclosure.md b/docs/src/tutorials/using_momentclosure.md
index 9341534..390e8bf 100644
--- a/docs/src/tutorials/using_momentclosure.md
+++ b/docs/src/tutorials/using_momentclosure.md
@@ -19,10 +19,11 @@ The terminology and the notation used throughout is consistent with the **Theory
using Catalyst
rn = @reaction_network begin
# including system-size parameter Ω
+ @parameters c₁ c₂ c₃ c₄ Ω
(c₁/Ω^2), 2X + Y → 3X
(c₂), X → Y
(c₃*Ω, c₄), 0 ↔ X
-end c₁ c₂ c₃ c₄ Ω
+end
```
The returned `rn` is an instance of [`ModelingToolkit.ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem). The net stoichiometry matrix and an array of the corresponding propensities, if needed, can be extracted directly from the model using [`Catalyst.netstoichmat`](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.netstoichmat) and MomentClosure function [`propensities`](@ref) respectively.
@@ -37,7 +38,7 @@ We can now obtain the moment equations. The system follows the law of mass actio
Let's start with the raw moment equations which we choose to generate up to second order ($m=2$):
```julia
using MomentClosure
-raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false)
+raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)
```
Note that we have set `combinatoric_ratelaw=false` in order to ignore the factorial scaling factors which [Catalyst adds to mass-action reactions](https://catalyst.sciml.ai/stable/tutorials/models/#Reaction-rate-laws-used-in-simulations). The function [`generate_raw_moment_eqs`](@ref) returns an instance of [`RawMomentEquations`](@ref MomentClosure.RawMomentEquations) that contains a [`ModelingToolkit.ODESystem`](https://mtk.sciml.ai/stable/systems/ODESystem/) composed of all the moment equations (accessed by `raw_eqs.odes`).
@@ -75,7 +76,7 @@ Coming back to the generated moment equations, we observe that they depend on hi
The corresponding central moment equations can also be easily generated:
```julia
-central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false)
+central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false)
```
Note that in case of non-polynomial propensity functions the [Taylor expansion order $q$](@ref central_moment_eqs) must also be specified, see the [P53 system example](P53_system_example.md) for more details. Luckily, the Brusselator contains only mass-action reactions and hence $q$ is automatically determined by the highest order (polynomial) propensity. The function [`generate_central_moment_eqs`](@ref) returns an instance of [`CentralMomentEquations`](@ref MomentClosure.CentralMomentEquations). As before, we can visualise the central moment equations:
```julia
@@ -189,7 +190,7 @@ using OrdinaryDiffEq
sol = solve(oprob, Tsit5(), saveat=0.1)
using Plots
-plot(sol, vars=(0, [1,2]), lw=2)
+plot(sol, idxs=[1,2], lw=2)
```
@@ -213,7 +214,7 @@ dprob = DiscreteProblem(jsys, u₀, tspan, p) # same parameters as defined earli
jprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false, false))
# define an EnsembleProblem to simulate multiple trajectories
-ensembleprob = EnsembleProblem(jprob)
+ensembleprob = EnsembleProblem(jprob)
# simulate 10000 SSA trajectories
@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.1, trajectories=10000)