diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..ff6499d --- /dev/null +++ b/.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" \ No newline at end of file diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 01e71aa..245daf2 100644 --- a/.github/workflows/Documentation.yml +++ b/.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 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0308eb3..7544a1b 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/Project.toml b/Project.toml index f44b30d..0759ce6 100644 --- a/Project.toml +++ b/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" 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") ``` ![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 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) ``` ![P53-Mdm2 distribution](../assets/p53-Mdm2_distribution.svg) @@ -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) ``` ![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") 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) ``` ![Brusselator issue 1](../assets/brusselator_issue_1.svg) @@ -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) ``` ![Brusselator issue 2](../assets/brusselator_issue_2.svg) @@ -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) ``` ![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. 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) ``` ![Time-dependent Brusselator normal](../assets/Brusselator_time-dependent_normal.svg) 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) ``` ![Brusselator means 1](../assets/brusselator_means_1.svg) @@ -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) diff --git a/src/MomentClosure.jl b/src/MomentClosure.jl index 448d118..1da2453 100644 --- a/src/MomentClosure.jl +++ b/src/MomentClosure.jl @@ -7,11 +7,12 @@ using SciMLBase, SciMLBase.EnsembleAnalysis using Random using Distributions: Geometric -using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, gradient, setmetadata, scalarize +using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, + gradient, setmetadata, scalarize using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint -using SymbolicUtils: Symbolic, Term, Real, expand, simplify, operation, - arguments, @rule, @acrule, isnotflat, flatten_term, - istree, FnType, Mul, Add, Pow, Div +using SymbolicUtils: BasicSymbolic, Real, Term, FnType, expand, simplify, + operation, arguments, @rule, @acrule, isnotflat, flatten_term, + istree, isterm, ismul, isadd, ispow, isdiv using DataStructures: OrderedDict using TupleTools: sort diff --git a/src/central_moment_equations.jl b/src/central_moment_equations.jl index 57f2809..7567445 100644 --- a/src/central_moment_equations.jl +++ b/src/central_moment_equations.jl @@ -17,7 +17,7 @@ end """ generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) Given a [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) return the [`CentralMomentEquations`](@ref) of the system generated up to `m_order`. @@ -32,7 +32,7 @@ Notes: - if `langevin=true`, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). *Note* that this field is irrelevant using `ReactionSystemMod` as then the @@ -42,13 +42,12 @@ Notes: accessible with [`Catalyst.speciesmap`](https://catalyst.sciml.ai/stable/api/catalyst_api/#Catalyst.speciesmap). """ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) N = numspecies(rn) # no. of molecular species in the network R = numreactions(rn) # no. of reactions in the network iv = get_iv(rn) # independent variable (usually time) - a = propensities(rn; combinatoric_ratelaw) # propensity functions of all reactions in the network - a = div_to_pow.(a) # convert Div to Pow (more stable at the moment than dealing with symbolic fractions) + a = propensities(rn; combinatoric_ratelaws) # propensity functions of all reactions in the network S = get_stoichiometry(rn, smap) # net stoichiometric matrix if langevin @@ -117,9 +116,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: suma = 0 for j in iter_all suma = suma + Da[r][j]*M[j]*1//fact(j) - #suma = simplify(suma) # here // gives fractions (otherwise it's a float number) - # only issue is that it does not latexify properly end for i in 1:N @@ -128,7 +125,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: else du[i] = S[i, r]*suma + du[i] end - du[i] = expand(du[i]) + du[i] = expand_expr(du[i]) end end @@ -158,7 +155,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: dM[i] -= i[j]*du[j]*M[i.-iter_1[j]] end end - dM[i] = expand(dM[i]) + dM[i] = expand_expr(dM[i]) end D = Differential(iv) diff --git a/src/closure_methods/closure.jl b/src/closure_methods/closure.jl index ceefc86..4d37ca4 100644 --- a/src/closure_methods/closure.jl +++ b/src/closure_methods/closure.jl @@ -1,7 +1,6 @@ function close_eqs(sys::MomentEquations, closure_exp::OrderedDict, closure::OrderedDict, polynorm::Bool) - # TODO: improve performance closed_eqs = Equation[] for eq in sys.odes.eqs @@ -10,7 +9,7 @@ function close_eqs(sys::MomentEquations, closure_exp::OrderedDict, if polynorm # depending on the functional form closed_rhs = expand(closed_rhs) else - closed_rhs = expand_expr(closed_rhs) + closed_rhs = expand_mod(closed_rhs) end push!(closed_eqs, Equation(eq.lhs, closed_rhs)) end diff --git a/src/closure_methods/conditional_derivative_matching.jl b/src/closure_methods/conditional_derivative_matching.jl index 6d664bb..09f3b9a 100644 --- a/src/closure_methods/conditional_derivative_matching.jl +++ b/src/closure_methods/conditional_derivative_matching.jl @@ -163,11 +163,10 @@ function conditional_derivative_matching(sys::MomentEquations, raw_to_central_exp = raw_to_central_moments(N, sys.q_order, μ_M_exp, bernoulli=true) raw_to_central = raw_to_central_moments(N, sys.q_order, μ_M, bernoulli=true) for i in sys.iter_q - closure_exp[sys.M[i]] = simplify(raw_to_central_exp[i]) - closure[sys.M[i]] = simplify(raw_to_central[i]) + closure_exp[sys.M[i]] = expand_mod(raw_to_central_exp[i]) + closure[sys.M[i]] = expand_mod(raw_to_central[i]) end - else closure_exp = closure_μ_exp closure = closure_μ diff --git a/src/closure_methods/derivative_matching.jl b/src/closure_methods/derivative_matching.jl index 15d6152..0f1c460 100644 --- a/src/closure_methods/derivative_matching.jl +++ b/src/closure_methods/derivative_matching.jl @@ -5,7 +5,7 @@ function multi_binomial(a, b) if length(a) != length(b) throw(ErrorException("vector lengths are different for some reason")) else - C = prod([binomial(a[i], b[i]) for i in 1:length(a)]) + C = prod([binomial(a[i], b[i]) for i in eachindex(a)]) return C end end diff --git a/src/closure_methods/gamma_closure.jl b/src/closure_methods/gamma_closure.jl index b2d9d09..545fc44 100644 --- a/src/closure_methods/gamma_closure.jl +++ b/src/closure_methods/gamma_closure.jl @@ -185,12 +185,12 @@ function gamma_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) for i in sys.iter_q closure_exp[sys.μ[i]] = substitute(μ[i], M_to_μ) expr = substitute(closure[sys.μ[i]], M_to_μ) - closure[sys.μ[i]] = simplify(expr) + closure[sys.μ[i]] = simplify(expr, simplify_fractions=false) end end if isempty(binary_vars) - return close_eqs(sys, closure_exp, closure, true) + return close_eqs(sys, closure_exp, closure, false) else # have to perform the bernoulli simplify in the end for gamma closure # as otherwise some symbolic terms do not cancel out and break it @@ -204,7 +204,7 @@ function gamma_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) closed_rhs = substitute(eq.rhs, closure_exp) closed_rhs = expand(closed_rhs) closed_rhs = substitute(closed_rhs, iter_sub) - closed_rhs = simplify(closed_rhs) + closed_rhs = simplify(closed_rhs, simplify_fractions=false) push!(closed_eqs, Equation(eq.lhs, closed_rhs)) diff --git a/src/closure_methods/linear_mapping_approximation.jl b/src/closure_methods/linear_mapping_approximation.jl index 0a5859c..43d32a3 100644 --- a/src/closure_methods/linear_mapping_approximation.jl +++ b/src/closure_methods/linear_mapping_approximation.jl @@ -1,6 +1,6 @@ """ linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0; - combinatoric_ratelaw = true) where T <: ReactionSystem + combinatoric_ratelaws = true) where T <: ReactionSystem Given a *nonlinear* [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) and an equivalent *linear* `ReactionSystem`, perform the Linear Mapping Approximation (LMA) @@ -21,12 +21,12 @@ Notes: - By default the moment equations will be generated up to the order determined by the degree of nonlinearity of the nonlinear system's reactions. However, if higher order moment information is required, the optional `m_order` argument may be provided to increase the expansion order manually. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). """ function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0; - combinatoric_ratelaw = true) where T <: ReactionSystem + combinatoric_ratelaws = true) where T <: ReactionSystem # check that all necessary information is provided and that LMA is applicable on the given nonlinear system @@ -47,11 +47,11 @@ function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars # apply LMA: substitute reaction parameters in the linear network moment equations to reflect # the nonlinear interactions according to the LMA methodology - term_factors, term_powers, poly_order = polynomial_propensities(propensities(rn_nonlinear; combinatoric_ratelaw)[nonlinear_rs_inds], + term_factors, term_powers, poly_order = polynomial_propensities(propensities(rn_nonlinear; combinatoric_ratelaws)[nonlinear_rs_inds], get_iv(rn_nonlinear), speciesmap(rn_nonlinear)) order = iszero(m_order) ? poly_order : max(poly_order, m_order) - try sys = generate_raw_moment_eqs(rn_linear, order; combinatoric_ratelaw, smap=speciesmap(rn_nonlinear)) + try sys = generate_raw_moment_eqs(rn_linear, order; combinatoric_ratelaws, smap=speciesmap(rn_nonlinear)) catch e error("LMA cannot handle reactions with non-polynomial rates\n $e") end diff --git a/src/closure_methods/log_normal_closure.jl b/src/closure_methods/log_normal_closure.jl index 0174eb3..e1ee438 100644 --- a/src/closure_methods/log_normal_closure.jl +++ b/src/closure_methods/log_normal_closure.jl @@ -1,7 +1,7 @@ function log_normal_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) closure = OrderedDict() - closure_exp = OrderedDict() # here it does not play a role + closure_exp = OrderedDict() N = sys.N if isempty(binary_vars) @@ -76,9 +76,9 @@ function log_normal_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[ central_to_raw = central_to_raw_moments(N, sys.q_order) closure_M = OrderedDict() for i in sys.iter_q - closure_exp[M[i]] = raw_to_central[i] + closure_exp[M[i]] = expand_mod(raw_to_central[i]) closure_M[M[i]] = closure[μ_symbolic[i]]-(central_to_raw[i]-M[i]) - closure_M[M[i]] = simplify(closure_M[M[i]]) + closure_M[M[i]] = expand_mod(closure_M[M[i]]) end closure = closure_M else diff --git a/src/moment_convert.jl b/src/moment_convert.jl index 6016606..dc4035d 100644 --- a/src/moment_convert.jl +++ b/src/moment_convert.jl @@ -183,7 +183,7 @@ function raw_to_central_moments(N::Int, order::Int, μ=nothing; bernoulli=false) end suma += term end - raw_to_central[i] = simplify(suma) + raw_to_central[i] = simplify(suma, simplify_fractions=false) end raw_to_central @@ -217,7 +217,7 @@ function central_to_raw_moments(N::Int, order::Int) end suma += term end - central_to_raw[i] = simplify(suma) + central_to_raw[i] = simplify(suma, simplify_fractions=false) end central_to_raw diff --git a/src/raw_moment_equations.jl b/src/raw_moment_equations.jl index 8390c19..2742dfa 100644 --- a/src/raw_moment_equations.jl +++ b/src/raw_moment_equations.jl @@ -1,6 +1,6 @@ """ generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) Given a [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) return the [`RawMomentEquations`](@ref) of the system generated up to `m_order`. @@ -13,7 +13,7 @@ Notes: - if `langevin=true`, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). *Note* that this field is irrelevant using `ReactionSystemMod` as then the @@ -23,13 +23,12 @@ Notes: accessible with [`Catalyst.speciesmap`](https://catalyst.sciml.ai/stable/api/catalyst_api/#Catalyst.speciesmap). """ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) iv = get_iv(rn) N = numspecies(rn) S = get_stoichiometry(rn, smap) - a = propensities(rn; combinatoric_ratelaw) - a = div_to_pow.(a) # NOTE: more stable at the moment than dealing with symbolic fractions + a = propensities(rn; combinatoric_ratelaws) if langevin drift = S*a @@ -72,7 +71,7 @@ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; dμ[i] += factor_j * suma end end - dμ[i] = expand(dμ[i]) + dμ[i] = expand_expr(dμ[i]) end D = Differential(iv) diff --git a/src/reaction_systems.jl b/src/reaction_systems.jl index d7fc320..11d2de2 100644 --- a/src/reaction_systems.jl +++ b/src/reaction_systems.jl @@ -1,15 +1,15 @@ """ - propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaw=true) + propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaws=true) Return a vector of propensity functions of all reactions in the given [`ReactionSystem`] (https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem). Notes: -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). """ -function propensities(rn::ReactionSystem; combinatoric_ratelaw=true) - simplify.(jumpratelaw.(reactions(rn); combinatoric_ratelaw)) +function propensities(rn::ReactionSystem; combinatoric_ratelaws::Bool=true) + simplify.(jumpratelaw.(reactions(rn); combinatoric_ratelaw=combinatoric_ratelaws)) end """ diff --git a/src/stochastic_stoichiometry.jl b/src/stochastic_stoichiometry.jl index 8e7538e..505b800 100644 --- a/src/stochastic_stoichiometry.jl +++ b/src/stochastic_stoichiometry.jl @@ -4,21 +4,31 @@ # TODO: follow changes in Catalyst and extend the compatibility accordingly # TODO: extend for other univariate discrete distributions -expected_coeff(x, n::Int) = x^n # covers Number/Sym/Pow -expected_coeff(x::Mul, n::Int) = prod(expected_coeff(a, n) for a in arguments(x)) +expected_coeff_mul(x) = prod(expected_coeff(a) for a in arguments(x)) +expected_coeff_pow(x) = expected_coeff(arguments(x)...) +expected_coeff_term(x) = expected_coeff(x, 1) -function expected_coeff(x::Div, n::Int) +function expected_coeff(x) + if ismul(x) + expected_coeff_mul(x) + elseif ispow(x) + expected_coeff_pow(x) + elseif isterm(x) + expected_coeff_term(x) + else + x + end +end + +expected_coeff_add(x, n::Int) = sum( expected_coeff(a) for a in arguments(expand(x^n)) ) +expected_coeff_mul(x, n::Int) = prod(expected_coeff(a, n) for a in arguments(x)) + +function expected_coeff_div(x, n::Int) num, denom = arguments(x) (expected_coeff(num, n) / expected_coeff(denom, n)) end -expected_coeff(x::Add, n::Int) = sum( expected_coeff(a) for a in arguments(expand(x^n)) ) -expected_coeff(x::Mul) = prod(expected_coeff(a) for a in arguments(x)) -expected_coeff(x::Pow) = expected_coeff(arguments(x)...) -expected_coeff(x::Term) = expected_coeff(x, 1) -expected_coeff(x) = x - -function expected_coeff(x::Term, n::Int) +function expected_coeff_term(x, n::Int) if isequal(operation(x), rand) args = arguments(x) isone(length(args)) || error("Unexpected arguments in $x") @@ -28,6 +38,21 @@ function expected_coeff(x::Term, n::Int) end end +function expected_coeff(x, n::Int) + if ismul(x) + expected_coeff_mul(x, n) + elseif isdiv(x) + expected_coeff_div(x, n) + elseif isadd(x) + expected_coeff_add(x, n) + elseif isterm(x) + expected_coeff_term(x, n) + else + x^n + end +end + + function resolve_moment(d, args, n) @assert isequal(d, Geometric) "Can only compute moments for the geometric distribution." geometric_raw_moment(args[1], n) diff --git a/src/symbolic.jl b/src/symbolic.jl index 83c9f15..270089c 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -19,14 +19,13 @@ end # Trim a string of form "(a, b, c, d, ...)" to "abcd..." trim_key(expr) = filter(x -> !(isspace(x) || x == ')' || x== '(' || x==','), string(expr)) -# Expand a symbolic expression (no binomial expansion) -expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) # apply distribution law -expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) # distributes terms until no longer possible -flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # -flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # -expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible +# Symbolic rules to expand expressions avoiding binomial expansion and fraction simplification +expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) +expand_expr = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) +expand_div = Fixpoint(Prewalk(PassThrough(@rule( +(~~xs) / ~a => sum(map(x -> x / ~a, ~~xs) ))))) +expand_mod = Fixpoint(Chain([expand_div, expand_expr])) -function define_μ(iter::AbstractVector, iv::Union{Sym,Num}) +function define_μ(iter::AbstractVector, iv::BasicSymbolic) indices = map(trim_key, iter) @@ -47,15 +46,15 @@ function define_μ(iter::AbstractVector, iv::Union{Sym,Num}) end -define_μ(N::Int, order::Int, iv::Union{Sym,Num}) = define_μ(construct_iter_all(N, order), iv) +define_μ(N::Int, order::Int, iv::BasicSymbolic) = define_μ(construct_iter_all(N, order), iv) -function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) +function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t return define_μ(iter, value(t)) end -function define_M(iter::AbstractVector, iv::Union{Sym,Num}) +function define_M(iter::AbstractVector, iv::BasicSymbolic) indices = map(trim_key, iter) @@ -78,15 +77,15 @@ function define_M(iter::AbstractVector, iv::Union{Sym,Num}) end -define_M(N::Int, order::Int, iv::Union{Sym,Num}) = define_M(construct_iter_all(N, order), iv) +define_M(N::Int, order::Int, iv::BasicSymbolic) = define_M(construct_iter_all(N, order), iv) -function define_M(N::Int, order::Int, iter = construct_iter_all(N, order)) +function define_M(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t return define_M(iter, value(t)) end function extract_variables(eqs::Array{Equation, 1}, μ, M=[]) - + vars = vcat(values(μ)..., values(M)...) # extract variables from rhs of each equation eq_vars = unique(vcat(get_variables.(eqs)...)) @@ -101,15 +100,6 @@ function extract_variables(eqs::Array{Equation, 1}, μ, M=[]) end -# a simple hack to convert SymbolicUtils Div (fractions) to Pow (multiplication) -# more stable in symbolic manipulations that we perform -function div_to_pow(expr::Div) - args = arguments(expr) - args[1] * args[2]^-1 -end - -div_to_pow(expr) = expr - ## Set of functions to deconstruct polynomial propensities ## #= @@ -133,21 +123,20 @@ isvar(x, vars) = any(isequal(x), vars) isconstant(expr, vars, iv) = !istree(expr) || (!isvar(expr, vars) && all(arg -> isconstant(arg, vars, iv), arguments(expr))) -function split_factor(expr::Pow, iv, smap, vars) +function split_factor_pow(expr, iv, vars) base, exp = arguments(expr) - #(exp isa Int && exp >= 0) || error("Unexpected exponent: $expr") exp isa Int || error("Unexpected exponent: $expr") - factor, powers = split_factor(base, iv, smap, vars) + factor, powers = split_factor(base, iv, vars) factor ^ exp, powers .* exp end -function split_factor(expr::Mul, iv, smap, vars) - powers = zeros(Int, length(smap)) +function split_factor_mul(expr, iv, vars) + powers = zeros(Int, length(vars)) factor = 1 for arg in arguments(expr) - factor_arg, power_arg = split_factor(arg, iv, smap, vars) + factor_arg, power_arg = split_factor(arg, iv, vars) factor *= factor_arg powers .+= power_arg end @@ -155,17 +144,24 @@ function split_factor(expr::Mul, iv, smap, vars) factor, powers end -function split_factor(expr::Div, iv, smap, vars) +function split_factor_div(expr, iv, vars) num, denom = arguments(expr) - isconstant(denom, vars, iv) || error("The denominator $denom is not constant.") - - factor, powers = split_factor(num, iv, smap, vars) - + isconstant(denom, vars, iv) || error("The denominator $denom in propensity $expr is not constant.") + + factor, powers = split_factor(num, iv, vars) + factor / denom, powers end -function split_factor(expr, iv, smap, vars) - if isconstant(expr, vars, iv) +function split_factor(expr, iv, vars) + + if ispow(expr) + split_factor_pow(expr, iv, vars) + elseif ismul(expr) + split_factor_mul(expr, iv, vars) + elseif isdiv(expr) + split_factor_div(expr, iv, vars) + elseif isconstant(expr, vars, iv) expr, zeros(length(vars)) elseif isvar(expr, vars) 1, map(isequal(expr), vars) @@ -174,20 +170,21 @@ function split_factor(expr, iv, smap, vars) end end -function polynomial_propensity(expr::Div, iv, smap, vars) + +function polynomial_propensity_div(expr, iv, vars) num, denom = arguments(expr) - isconstant(denom, vars, iv) || error("The denominator $denom is not constant.") - factors, powers = polynomial_propensity(num, iv, smap, vars) - factors ./ denom, powers + isconstant(denom, vars, iv) || error("The denominator $denom in propensity $expr is not constant.") + factors, powers = polynomial_propensity(num, iv, vars) + factors ./denom, powers end -function polynomial_propensity(expr::Add, iv, smap, vars) +function polynomial_propensity_add(expr, iv, vars) factors = [] powers = Vector{Int}[] - + for term in arguments(expr) factor_term, power_term = try - split_factor(term, iv, smap, vars) + split_factor(term, iv, vars) catch e error("Propensity function $term is non-polynomial? \n" * string(e)) end @@ -198,17 +195,18 @@ function polynomial_propensity(expr::Add, iv, smap, vars) factors, powers end -function polynomial_propensity(expr, iv, smap, vars) - factor, powers = try - split_factor(expr, iv, smap, vars) - catch e - error("Propensity function $expr is non-polynomial? \n" * string(e)) +function polynomial_propensity(expr, iv, vars) + if isdiv(expr) + polynomial_propensity_div(expr, iv, vars) + elseif isadd(expr) + polynomial_propensity_add(expr, iv, vars) + else + factor, powers = split_factor(expr, iv, vars) + [factor], [powers] end - - [ factor ], [ powers ] end -function polynomial_propensities(arr::AbstractArray, iv::Sym, smap::AbstractDict) +function polynomial_propensities(arr::AbstractArray, iv::BasicSymbolic, smap::AbstractDict) vars = [ x for (x,_) in Base.sort(collect(smap), by=x->x[2]) ] all_factors = Array{Vector}(undef, size(arr)) @@ -216,9 +214,9 @@ function polynomial_propensities(arr::AbstractArray, iv::Sym, smap::AbstractDict for (rind, expr) in enumerate(arr) expr = expand(expr) - all_factors[rind], all_powers[rind] = polynomial_propensity(expr, iv, smap, vars) + all_factors[rind], all_powers[rind] = polynomial_propensity(expr, iv, vars) end max_power = maximum(sum.(vcat(all_powers...))) all_factors, all_powers, max_power -end \ No newline at end of file +end diff --git a/src/utils.jl b/src/utils.jl index 1c611c9..99053fc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -307,6 +307,7 @@ function format_moment_eqs(eqs::MomentEquations) for i in 1:size(odes)[1] key = vars[i] eq = odes[i].rhs + eq = expand_mod(eq) expr = "d"*string(key)*"/dt = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") @@ -341,6 +342,7 @@ function format_closure(eqs::ClosedMomentEquations; format_all::Bool=false) for i in iter eq = closure[i] + eq = expand_mod(eq) expr = string(i)*" = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") diff --git a/test/bernoulli_closure_tests.jl b/test/bernoulli_closure_tests.jl index ae25387..ac75972 100644 --- a/test/bernoulli_closure_tests.jl +++ b/test/bernoulli_closure_tests.jl @@ -13,13 +13,13 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] μ = define_μ(2,4) M = define_M(2,4) -sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) closed_eqs = moment_closure(sys, "zero", binary_vars) @test length(closed_eqs.odes.states) == 4 && isequal(closed_eqs.closure[M[0,4]], 0) @@ -30,7 +30,7 @@ closed_eqs = moment_closure(sys, "normal", binary_vars) closed_eqs= moment_closure(sys, "log-normal", binary_vars) expr1 = expand(closed_eqs.closure[M[0,3]]) expr2 = M[0,2]^3*μ[0,1]^-3 + 3*M[0,2]^2*μ[0,1]^-1 -@test length(closed_eqs.odes.states) == 4 && isequal(expand(expr1), expr2) +@test length(closed_eqs.odes.states) == 4 && isequal(simplify(expr1), simplify(expr2)) closed_eqs = moment_closure(sys, "poisson", binary_vars) @test length(closed_eqs.odes.states) == 4 && isequal(closed_eqs.closure[M[0,3]], μ[0,1]) @@ -43,7 +43,7 @@ expr2 = 2*M[0,2]^2*μ[0,1]^-1 closed_eqs = moment_closure(sys, "derivative matching", binary_vars) expr1 = expand(closed_eqs.closure[sys.M[0,3]]) expr2 = M[0,2]^3*μ[0,1]^-3 + 3*M[0,2]^2*μ[0,1]^-1 -@test length(closed_eqs.odes.states) == 4 && isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 && isequal(simplify(expr1), simplify(expr2)) sys = generate_raw_moment_eqs(rn, 2) @@ -60,7 +60,7 @@ closed_eqs = moment_closure(sys, "poisson", binary_vars) @test length(closed_eqs.odes.states) == 4 closed_eqs = moment_closure(sys, "gamma", binary_vars) -@test isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 closed_eqs = moment_closure(sys, "derivative matching", binary_vars) -@test isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 \ No newline at end of file diff --git a/test/closure_tests.jl b/test/closure_tests.jl index adf8d73..0f014b6 100644 --- a/test/closure_tests.jl +++ b/test/closure_tests.jl @@ -1,6 +1,7 @@ using MomentClosure using MomentClosure: define_M, define_μ using Symbolics: value, expand +using SymbolicUtils: Fixpoint using Test using Catalyst @@ -8,24 +9,22 @@ rn = @reaction_network begin (c₁/Ω^2), 2X + Y → 3X (c₂), X → Y (Ω*c₃, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end -#@parameters c₁, c₂, c₃, c₄, Ω @syms c₁::Real c₂::Real c₃::Real c₄::Real Ω::Real μ = define_μ(2,4) M = define_M(2,4) # --- Test closures on central moment equations --- -sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[1].rhs closed_eqs = moment_closure(sys, "zero") @test closed_eqs.closure[M[2,2]] == 0 expr1 = closed_eqs.odes.eqs[1].rhs expr2 = c₃*Ω + M[1,1]*c₁*μ[1,0]*(Ω^-2) + M[1,1]*c₁*(Ω^-2)*(μ[1,0]- 1) + c₁*M[2,0]*μ[0,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*(Ω^-2)*(μ[1,0] - 1) - c₂*μ[1,0] - c₄*μ[1,0] -expr2 = value.(expand(expr2)) -@test isequal(expand(expr1), expand(expr2)) +@test isequal(simplify(expr1), simplify(expr2)) # check that deterministic_IC is working with central moments ic_values = Dict(deterministic_IC([2, 5], closed_eqs)) @@ -39,7 +38,7 @@ closed_eqs= moment_closure(sys, "log-normal") expr1 = closed_eqs.closure[M[1,2]] expr2 = μ[1,0]*μ[0,1]^2*(1.0+M[0,2]*μ[0,1]^-2)*(1.0 + M[1,1]*(μ[0,1]^-1)*(μ[1,0]^-1))^2 - M[0,2]*μ[1,0] - μ[1,0]*μ[0,1]^2 - 2*M[1,1]*μ[0,1] -@test isequal(expand(expr1), expand(expr2)) +@test isequal(Fixpoint(simplify)(expr1), Fixpoint(simplify)(expr2)) closed_eqs = moment_closure(sys, "poisson") @test isequal(closed_eqs.closure[M[3,0]], μ[1,0]) @@ -54,11 +53,11 @@ closed_eqs = moment_closure(sys, "derivative matching") expr1 = closed_eqs.closure[sys.M[0,4]] expr2 = μ[0,1]^4*(M[0,2]+μ[0,1]^2)^-6*(M[0,3]+μ[0,1]^3+3*M[0,2]*μ[0,1])^4 - μ[0,1]^4 - 6*M[0,2]*μ[0,1]^2 - 4*M[0,3]*μ[0,1] -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), simplify(expr2)) # --- Test closures on raw moment equations --- -sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) closed_eqs = moment_closure(sys, "zero") expr1 = closed_eqs.closure[μ[3,0]] expr2 = -2*μ[1,0]^3 + 3*μ[1,0]*μ[2,0] @@ -72,8 +71,7 @@ expr1 = closed_eqs.odes.eqs[5].rhs expr2 = c₂*μ[1,0] + 2*c₂*μ[1,1] + c₁*μ[0,1]*μ[2,0]*Ω^-2 - c₁*μ[1,1]*Ω^-2 - 4*c₁*μ[1,1]^2*Ω^-2 + 4*c₁*μ[0,1]*μ[1,1]*Ω^-2 - 2*c₁*μ[0,1]*μ[1,0]^2*Ω^-2 + 2*c₁*μ[0,2]*μ[1,0]*Ω^-2 -2*c₁*μ[0,2]*μ[2,0]*Ω^-2 + 2*c₁*μ[1,0]*μ[1,1]*Ω^-2 - 4*c₁*μ[1,0]*μ[0,1]^2*Ω^-2 + 4*c₁*μ[0,1]^2*μ[1,0]^2*Ω^-2 -expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expr2) +@test isequal(simplify(expr1), simplify(expr2)) # check that deterministic_IC is working with raw moments ic_values = Dict(deterministic_IC([2, 5], closed_eqs)) @@ -93,7 +91,7 @@ expr2 = μ[1,0] + 6 * μ[1,0]^4 + 3*μ[2,0]^2 + 4*μ[1,0]*μ[3,0] -12*μ[2,0]*μ closed_eqs = moment_closure(sys, "gamma") expr1 = closed_eqs.closure[μ[0,3]] expr2 = μ[0,1]^3 + 3*μ[0,1]*(μ[0,2]-μ[0,1]^2) + 2*μ[0,1]^-1*(μ[0,2]-μ[0,1]^2)^2 -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), simplify(expr2)) closed_eqs = moment_closure(sys, "derivative matching") expr1 = closed_eqs.closure[sys.μ[0,3]] diff --git a/test/conditional_closure_tests.jl b/test/conditional_closure_tests.jl index bca71c8..3134689 100644 --- a/test/conditional_closure_tests.jl +++ b/test/conditional_closure_tests.jl @@ -14,7 +14,7 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] @parameters k_on, k_off, k_p, γ_p @@ -39,10 +39,10 @@ expr2 = simplify(value.(expr2)) closed_eqs = moment_closure(sys, "conditional gaussian", binary_vars) expr1 = closed_eqs.closure[μ[1,4]] expr2 = 4*μ[1,1]*μ[1,3]*μ[1,0]^-1 + 3*μ[1,2]^2*μ[1,0]^-1 - 12*μ[1,2]*μ[1,1]^2*μ[1,0]^-2 + 6*μ[1,1]^4*μ[1,0]^-3 -@test isequal(expr1, expr2) +@test isequal(expr1, simplify(expr2)) expr1 = closed_eqs.closure[μ[1,3]] expr2 = 3*μ[1,2]*μ[1,1]*μ[1,0]^-1 - 2*μ[1,1]^3*μ[1,0]^-2 -@test isequal(expr1, expr2) +@test isequal(expr1, simplify(expr2)) closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) expr1 = closed_eqs.closure[μ[1,4]] @@ -70,8 +70,8 @@ expr1 = expand(closed_eqs.closure[M[0,5]]) expr2 = 10*M[0,2]*M[0,3] + 5*M[0,4]*μ[0,1] - 15*μ[0,1]*M[0,2]^2 @test isequal(expr1, expr2) -closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) +@time closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) expr1 = closed_eqs.closure[M[1,3]] expr2 = μ[1,0]*(M[1,1]+μ[0,1]*μ[1,0])^-3*(M[1,2]+M[0,2]*μ[1,0]+μ[1,0]*μ[0,1]^2+2*M[1,1]*μ[0,1])^3 - M[0,3]*μ[1,0] - 3*M[1,1]*μ[0,1]^2 - 3*M[1,2]*μ[0,1] - μ[1,0]*μ[0,1]^3- 3*M[0,2]*μ[0,1]*μ[1,0] -@test isequal(expand(expr1), expand(expr2)) +@test isequal(simplify(expr1), simplify(expr2)) \ No newline at end of file diff --git a/test/latexify_tests.jl b/test/latexify_tests.jl index 86687fa..002892f 100644 --- a/test/latexify_tests.jl +++ b/test/latexify_tests.jl @@ -13,10 +13,10 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) clean_eqs = bernoulli_moment_eqs(raw_eqs, binary_vars) closed_raw_eqs = moment_closure(raw_eqs, "conditional gaussian", binary_vars) @@ -37,11 +37,18 @@ expr = replace(raw"\begin{align*} \frac{d\mu_{1 0}}{dt} =& k_{on}", "\r\n"=>"\n") @test latexify(clean_eqs)[1:46] == expr +#= expr = replace(raw"\begin{align*} \mu_{1 2} =& \mu_{1 0}^{-1} \mu_{1 1}^{2} \\ \mu_{1 3} =& 3 \mu_{1 0}^{-1} \mu_{1 1} \mu_{1 2} - 2 \mu_{1 0}^{-2} \mu_{1 1}^{3} \end{align*} ", "\r\n"=>"\n") +=# +expr = replace(raw"\begin{align*} +\mu_{1 2} =& \frac{\mu_{1 1}^{2}}{\mu_{1 0}} \\ +\mu_{1 3} =& \frac{-2 \mu_{1 1}^{3}}{\mu_{1 0}^{2}} + \frac{3 \mu_{1 1} \mu_{1 2}}{\mu_{1 0}} +\end{align*} +", "\r\n"=>"\n") @test latexify(closed_raw_eqs, :closure) == expr @test_throws MethodError latexify(raw_eqs, :closure) \ No newline at end of file diff --git a/test/lma_tests.jl b/test/lma_tests.jl index 59c90a3..1aa27e8 100644 --- a/test/lma_tests.jl +++ b/test/lma_tests.jl @@ -14,7 +14,7 @@ rn_nonlinear = @reaction_network begin ρ_u, g → g + p ρ_b*(1-g), 0 ⇒ p 1, p → 0 -end σ_b σ_u ρ_b ρ_u +end rn_linear = @reaction_network begin σ_b_LMA, g → 0 @@ -22,10 +22,10 @@ rn_linear = @reaction_network begin ρ_u, g → g+p (ρ_b*(1-g)), 0 ⇒ p 1, p → 0 -end σ_b_LMA σ_u ρ_b ρ_u +end binary_vars = [speciesmap(rn_nonlinear)[g]] -LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) +LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false) μ = define_μ(2,3) expr1 = LMA_eqs.odes.eqs[1].rhs @@ -36,7 +36,7 @@ expr2 = ρ_b + ρ_u*μ[1,0] - ρ_b*μ[1,0] - μ[0,1] @test isequal(expr1, value.(expr2)) expr1 = LMA_eqs.odes.eqs[3].rhs expr2 = ρ_u*μ[1,0] + σ_u*μ[0,1] - μ[1,1] - σ_u*μ[1,1] - σ_b*μ[1,1]^2*μ[1,0]^-1 -@test isequal(expr1, value.(expr2)) +@test isequal(simplify(expr1), simplify(expr2)) # cooperativity cp=2 rn_nonlinear = @reaction_network begin @@ -45,10 +45,10 @@ rn_nonlinear = @reaction_network begin ρ_u, g → g + p ρ_b*(1-g), 0 ⇒ p 1, p → 0 -end σ_b σ_u ρ_b ρ_u +end binary_vars = [speciesmap(rn_nonlinear)[g]] -_, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) +_, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false) expr1 = effective_params[value(σ_b_LMA)] expr2 = (σ_b*μ[1,2] - σ_b*μ[1,1]) * μ[1,0]^-1 -@test isequal(expr1, value.(expr2)) +@test isequal(expr1, value(expr2)) \ No newline at end of file diff --git a/test/moment_equations_SDE_tests.jl b/test/moment_equations_SDE_tests.jl index 26f9c29..a1a85b2 100644 --- a/test/moment_equations_SDE_tests.jl +++ b/test/moment_equations_SDE_tests.jl @@ -24,15 +24,16 @@ cir_central_moments = generate_central_moment_eqs(cir_model, 2) ## Chemical Reaction Networks via Chemical Langevin Equation # unimolecular system schloegl = @reaction_network begin + @parameters k1 k2 1.0, 2X → 3X 1.0, 3X → 2X k1, ∅ → X k2, X → ∅ -end k1 k2 +end k1, k2 = parameters(schloegl) for order in 2:6 - schloegl_moments = generate_raw_moment_eqs(schloegl, order; langevin = true, combinatoric_ratelaw = false) + schloegl_moments = generate_raw_moment_eqs(schloegl, order; langevin = true, combinatoric_ratelaws = false) local μ = schloegl_moments.μ rhs(j) = j > 1 ? expand( j * ( k1*μ[(j-1,)] - k2*μ[(j,)] + (μ[(j+1,)] - μ[(j,)]) - (μ[(j+2,)] - 3*μ[(j+1,)] + 2*μ[(j,)]) ) + j*(j-1)/2 * ( k1*μ[(j-2,)] + k2*μ[(j-1,)] + (μ[(j,)] - μ[(j-1,)]) + (μ[(j+1,)] - 3*μ[(j,)] + 2*μ[(j-1,)]) ) ) : @@ -42,7 +43,7 @@ for order in 2:6 end μ = define_μ(1, 4); M = define_M(1, 4) -schloegl_moments = generate_central_moment_eqs(schloegl, 2, 4, langevin=true, combinatoric_ratelaw = false) +schloegl_moments = generate_central_moment_eqs(schloegl, 2, 4, langevin=true, combinatoric_ratelaws = false) expr = k1 + 4*μ[(1,)]^2 + 4*M[(2,)] - μ[(1,)]^3 - M[(3,)] - 3*μ[(1,)] - k2*μ[(1,)] - 3*M[(2,)]*μ[(1,)] @test isequal(schloegl_moments.odes.eqs[1].rhs, expr) expr = k1 + 9*M[(3,)] + k2*μ[(1,)] + μ[(1,)]^3 + 19*M[(2,)]*μ[(1,)] + μ[(1,)] - 2*μ[(1,)]^2 - @@ -51,15 +52,16 @@ expr = k1 + 9*M[(3,)] + k2*μ[(1,)] + μ[(1,)]^3 + 19*M[(2,)]*μ[(1,)] + μ[(1,) # check if things work for multiple species rn = @reaction_network begin + @parameters c1 c2 c3 c4 (c1), 2X + Y → 3X (c2), X → Y (c3, c4), 0 ↔ X -end c1 c2 c3 c4 +end c1, c2, c3, c4 = parameters(rn) for order in 2:10 - rn_moments = generate_raw_moment_eqs(rn, order; langevin = true, combinatoric_ratelaw = false) + rn_moments = generate_raw_moment_eqs(rn, order; langevin = true, combinatoric_ratelaws = false) local μ = rn_moments.μ rhs(i,j) = (i >= 1 ? i*( c1*(μ[(i+1,j+1)]-μ[(i,j+1)]) - (c2 + c4)*μ[(i,j)] + c3*μ[(i-1,j)] ) : 0) + (j >= 1 ? j*( -c1*(μ[(i+2,j)]-μ[(i+1,j)]) + c2*μ[(i+1,j-1)] ) : 0) + @@ -73,8 +75,8 @@ for order in 2:10 @test isequal(analytic_moment_eqs, rn_moments.odes.eqs) end -raw_eqs = generate_raw_moment_eqs(rn, 2, langevin = true, combinatoric_ratelaw=false) -central_eqs = generate_central_moment_eqs(rn, 2, langevin=true, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, langevin = true, combinatoric_ratelaws=false) +central_eqs = generate_central_moment_eqs(rn, 2, langevin=true, combinatoric_ratelaws=false) μ = define_μ(2, 4); M = define_μ(2, 4) central_to_raw = central_to_raw_moments(2, 4) subdict = Dict( μ[iter] => central_to_raw[iter] for iter in filter(x -> sum(x) > 1, raw_eqs.iter_all)) diff --git a/test/moment_equations_tests.jl b/test/moment_equations_tests.jl index 24b0612..cee2033 100644 --- a/test/moment_equations_tests.jl +++ b/test/moment_equations_tests.jl @@ -1,5 +1,5 @@ using MomentClosure -using MomentClosure: expected_coeff +using MomentClosure: expected_coeff, expand_div using Symbolics: value, simplify, expand using Distributions: Geometric using Catalyst @@ -17,8 +17,8 @@ b = value(b); d = value(d); p = value(p) @test isequal(expected_coeff(sin(b), 2), sin(b)^2) @register_symbolic Geometric(p) m = rand(Geometric(p)) -@test isequal(expected_coeff(m, 1), p^-1 - 1) -@test isequal(expected_coeff(m, 3), -1 + 7p^-1 - 12p^-2 + 6p^-3) +@test isequal(expand_div(expected_coeff(m, 1)), p^-1 - 1) +@test isequal(expand_div(expected_coeff(m, 3)), expand(-1 + 7/p - 12/p^2 + 6/p^3)) @register_symbolic Bernoulli(p) @test_throws AssertionError expected_coeff(rand(Bernoulli(p)), 1) @@ -30,9 +30,9 @@ rn = @reaction_network begin (c₁*Ω^-2), 2X + Y → 3X (c₂), X → Y (Ω*c₃, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end -sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) μ = sys.μ @@ -41,20 +41,21 @@ M = sys.M expr2 = c₂*μ[1,0] + c₁*M[1,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*Ω^-2 - c₁*M[2,1]*Ω^-2 - 2*c₁*M[1,1]*μ[1,0]*Ω^-2 - c₁*M[2,0]*μ[0,1]*Ω^-2 - c₁*μ[0,1]*Ω^-2*μ[1,0]^2 expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expr2) +@test isequal(simplify(expr1), expr2) + -sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) -@test isequal(expand(expr1), expr2) +@test isequal(simplify(expr1), expr2) -sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) μ = sys.μ expr1 = sys.odes.eqs[4].rhs expr2 = c₂*μ[2,0] + c₁*μ[1,1]*Ω^-2 - c₁*μ[3,1]*Ω^-2 -c₂*(μ[1,0] + μ[1,1]) + c₃*μ[0,1]*Ω - c₄*μ[1,1] - c₁*μ[1,2]*Ω^-2 + c₁*μ[2,2]*Ω^-2 expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expand(expr2)) +@test isequal(simplify(expr1), expr2) # corner case - a linear propensity with rate coefficient 1 rn = @reaction_network begin @@ -68,21 +69,21 @@ sys = generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (X^1.5), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*Y^3*X^1.5), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*sin(exp(X)+Y)), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*sin(exp(t))+Y^(3+1)), X ⇒ Y -end c₁ +end eqs = generate_raw_moment_eqs(rn, 2) @test eqs.q_order == 5 \ No newline at end of file diff --git a/test/reaction_system_tests.jl b/test/reaction_system_tests.jl index 97bd7cb..0c8fd2f 100644 --- a/test/reaction_system_tests.jl +++ b/test/reaction_system_tests.jl @@ -3,19 +3,20 @@ using MomentClosure using Catalyst @parameters t, c₁, c₂, c₃, c₄, Ω -@variables X(t), Y(t) +@species X(t), Y(t) S_mat = [ 1 -1 1 -1; -1 1 0 0] -# NOTE: only holds if combinatoric_ratelaw=false +# NOTE: only holds if combinatoric_ratelaws=false a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X] 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 smap = speciesmap(rn) @test isequal(get_stoichiometry(rn, smap), S_mat) -@test isequal(propensities(rn, combinatoric_ratelaw=false), a) \ No newline at end of file +@test isequal(propensities(rn, combinatoric_ratelaws=false), a)