update docs

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")
 ```
 ![LMA feedback loop mean protein number](../assets/LMA_feedback_loop_mean_protein_number.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

docs/src/tutorials/P53_system_example.md

@@ -40,21 +40,22 @@ 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:
 ```julia
 # → 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)
 ```
 ![P53-Mdm2 distribution](../assets/p53-Mdm2_distribution.svg)
 
@@ -120,15 +121,15 @@ 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)
 
         u₀map = deterministic_IC(u₀, closed_eqs)
         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)
 ```
 ![P53-Mdm2 normal means 2nd order expansion](../assets/p53-Mdm2_normal_2nd_order.svg)
 
 Log-normal closure:
 ```julia
-plot(plts[2], size=(750, 450), leftmargin=2mm)
+plot(plts[2], size=(750, 450), leftmargin=2PlotMeasures.mm)
 ```
 ![P53-Mdm2 log-normal means 2nd order expansion](../assets/p53-Mdm2_log-normal_2nd_order.svg)
 
@@ -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)
 ```
 ![P53-Mdm2 gamma means 2nd order expansion](../assets/p53-Mdm2_gamma_2nd_order.svg)
 
@@ -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")

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)
 ```
 ![Brusselator issue 1](../assets/brusselator_issue_1.svg)
 
@@ -40,74 +41,74 @@ 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)
 ```
 ![Brusselator issue 2](../assets/brusselator_issue_2.svg)
 
 We observe sustained oscillatory behaviour instead of the expected damped oscillations. This result is unphysical: single SSA trajectories (that may display sustained oscillations) get dephased over time and hence the ensemble average should always show damped or overdamped oscillations [1].
 
 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)
 ```
 ![Brusselator issue 3](../assets/brusselator_issue_3.svg)
 
 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)
 ```
 ![Brusselator issue 4](../assets/brusselator_issue_4.svg)
 
 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)
 ```
 ![Brusselator issue 5](../assets/brusselator_issue_5.svg)
 
 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)
 ```
 ![Brusselator issue 6](../assets/brusselator_issue_6.svg)
 
 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.

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))