Merge pull request #48 from ChrisRackauckas/patch-1

Check Catalyst 13

Author: augustinas1

.github/dependabot.yml

@@ -0,0 +1,7 @@
+# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates
+version: 2
+updates:
+  - package-ecosystem: "github-actions"
+    directory: "/" # Location of package manifests
+    schedule:
+      interval: "weekly"

.github/workflows/Documentation.yml

@@ -11,7 +11,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
@@ -20,5 +20,5 @@ jobs:
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
-          #DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
         run: julia --project=docs/ docs/make.jl

.github/workflows/ci.yml

@@ -16,12 +16,12 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
         with:
@@ -34,6 +34,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v3
         with:
           file: lcov.info

Project.toml

@@ -1,7 +1,7 @@
 name = "MomentClosure"
 uuid = "01a1b25a-ecf0-48c5-ae58-55bfd5393600"
 authors = ["Augustinas Sukys"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
@@ -19,17 +19,17 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
 [compat]
-Catalyst = "11, 12"
+Catalyst = "13"
 Combinatorics = "1"
 Cumulants = "1"
 DataStructures = "0.18"
-Distributions = "0.23, 0.24, 0.25"
+Distributions = "0.25"
 DocStringExtensions = "0.8, 0.9"
-Latexify = "0.15.14"
-ModelingToolkit = "6.5, 7, 8"
+Latexify = "0.15.14, 0.16"
+ModelingToolkit = "8"
 SciMLBase = "1.8.3"
-SymbolicUtils = "0.19"
-Symbolics = "3, 4"
+SymbolicUtils = "1"
+Symbolics = "5"
 TupleTools = "1.2.0"
 julia = "1.6"
 

docs/src/assets/derivative_matching_cumulant.svg

@@ -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/LMA_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")