Merge remote-tracking branch 'origin/master'

vyudu · vyudu · commit b8d020158f65 · 2024-06-05T09:54:44.000-04:00
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -14,7 +14,7 @@ jobs:
         group:
           - Core
         version:
-          - '1.10.2'
+          - '1'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.10.2'
+          version: '1'
       - name: Install dependencies
         run: julia --project=docs/ -e 'ENV["JULIA_PKG_SERVER"] = ""; using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/Project.toml b/Project.toml
@@ -42,6 +42,7 @@ DataStructures = "0.18"
 Combinatorics = "1.0.2"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
+DynamicPolynomials = "0.5"
 DynamicQuantities = "0.13.2"
 Graphs = "1.4"
 HomotopyContinuation = "2.9"
@@ -59,7 +60,7 @@ StructuralIdentifiability = "0.5.1"
 SymbolicUtils = "1.0.3"
 Symbolics = "5.27"
 Unitful = "1.12.4"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -37,7 +37,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 BenchmarkTools = "1.5"
-BifurcationKit = "0.3"
+BifurcationKit = "0.3.4"
 CairoMakie = "0.12"
 Catalyst = "13"
 DataFrames = "1.6"
diff --git a/docs/pages.jl b/docs/pages.jl
@@ -9,10 +9,10 @@ pages = Any[
         "model_creation/dsl_basics.md",
         "model_creation/dsl_advanced.md",
         #"model_creation/programmatic_CRN_construction.md",
-        #"model_creation/compositional_modeling.md",
+        "model_creation/compositional_modeling.md",
         #"model_creation/constraint_equations.md",
         # Events.
-        #"model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions.
+        "model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions.
         "model_creation/model_file_loading_and_export.md",# Distributed parameters, rates, and initial conditions.
         # Loading and writing models to files.
         "model_creation/model_visualisation.md",
@@ -21,7 +21,7 @@ pages = Any[
         "Model creation examples" => Any[
             "model_creation/examples/basic_CRN_library.md",
             "model_creation/examples/programmatic_generative_linear_pathway.md",
-            #"model_creation/examples/hodgkin_huxley_equation.md",
+            "model_creation/examples/hodgkin_huxley_equation.md",
             #"model_creation/examples/smoluchowski_coagulation_equation.md"
         ]
     ],
@@ -66,6 +66,6 @@ pages = Any[
     #     # Contributor's guide.
     #     # Repository structure.
     # ],
-    #"FAQs" => "faqs.md",
-    #"API" => "api.md"
+    "FAQs" => "faqs.md",
+    "API" => "api.md"
 ]
diff --git a/docs/src/api.md b/docs/src/api.md
@@ -35,37 +35,41 @@ corresponding chemical reaction ODE models, chemical Langevin equation SDE
 models, and stochastic chemical kinetics jump process models.
 
 ```@example ex1
-using Catalyst, DifferentialEquations, Plots
+using Catalyst, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Plots
 t = default_t()
 @parameters β γ
 @species S(t) I(t) R(t)
 
 rxs = [Reaction(β, [S,I], [I], [1,1], [2])
        Reaction(γ, [I], [R])]
 @named rs = ReactionSystem(rxs, t)
+rs = complete(rs)
 
 u₀map    = [S => 999.0, I => 1.0, R => 0.0]
 parammap = [β => 1/10000, γ => 0.01]
 tspan    = (0.0, 250.0)
 
 # solve as ODEs
 odesys = convert(ODESystem, rs)
+odesys = complete(odesys)
 oprob = ODEProblem(odesys, u₀map, tspan, parammap)
 sol = solve(oprob, Tsit5())
 p1 = plot(sol, title = "ODE")
 
 # solve as SDEs
 sdesys = convert(SDESystem, rs)
+sdesys = complete(sdesys)
 sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
-sol = solve(sprob, EM(), dt=.01)
+sol = solve(sprob, EM(), dt=.01, saveat = 2.0)
 p2 = plot(sol, title = "SDE")
 
 # solve as jump process
 jumpsys = convert(JumpSystem, rs)
+jumpsys = complete(jumpsys)
 u₀map    = [S => 999, I => 1, R => 0]
 dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
-jprob = JumpProblem(jumpsys, dprob, Direct())
-sol = solve(jprob, SSAStepper())
+jprob = JumpProblem(jumpsys, dprob, Direct(); save_positions = (false,false))
+sol = solve(jprob, SSAStepper(), saveat = 2.0)
 p3 = plot(sol, title = "jump")
 
 plot(p1, p2, p3; layout = (3,1))
diff --git a/docs/src/faqs.md b/docs/src/faqs.md
@@ -5,19 +5,17 @@ One can directly use symbolic variables to index into SciML solution objects.
 Moreover, observables can also be evaluated in this way. For example,
 consider the system
 ```@example faq1
-using Catalyst, DifferentialEquations, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 rn = @reaction_network ABtoC begin
   (k₊,k₋), A + B <--> C
 end
-
-# initial condition and parameter values
-setdefaults!(rn, [:A => 1.0, :B => 2.0, :C => 0.0, :k₊ => 1.0, :k₋ => 1.0])
 nothing    # hide
 ```
 Let's convert it to a system of ODEs, using the conservation laws of the system
 to eliminate two of the species:
 ```@example faq1
 osys = convert(ODESystem, rn; remove_conserved = true)
+osys = complete(osys)
 ```
 Notice the resulting ODE system has just one ODE, while algebraic observables
 have been added for the two removed species (in terms of the conservation law
@@ -28,7 +26,9 @@ observed(osys)
 Let's solve the system and see how to index the solution using our symbolic
 variables
 ```@example faq1
-oprob = ODEProblem(osys, [], (0.0, 10.0), [])
+u0 = [osys.A => 1.0, osys.B => 2.0, osys.C => 0.0]
+ps = [osys.k₊ => 1.0, osys.k₋ => 1.0]
+oprob = ODEProblem(osys, u0, (0.0, 10.0), ps)
 sol = solve(oprob, Tsit5())
 ```
 Suppose we want to plot just species `C`, without having to know its integer
@@ -44,8 +44,8 @@ sol[C]
 ```
 To evaluate `C` at specific times and plot it we can just do
 ```@example faq1
-t = range(0.0, 10.0, length=101)
-plot(t, sol(t, idxs = C), label = "C(t)", xlabel = "t")
+t = range(0.0, 10.0, length = 101)
+plot(sol(t, idxs = C), label = "C(t)", xlabel = "t")
 ```
 If we want to get multiple variables we can just do
 ```@example faq1
@@ -65,7 +65,7 @@ constant, giving `k*X^2/2` instead of `k*X^2` for ODEs and `k*X*(X-1)/2` instead
 of `k*X*(X-1)` for jumps. This can be disabled when directly `convert`ing a
 [`ReactionSystem`](@ref). If `rn` is a generated [`ReactionSystem`](@ref), we can
 do
-```julia
+```@example faq1
 osys = convert(ODESystem, rn; combinatoric_ratelaws=false)
 ```
 Disabling these rescalings should work for all conversions of `ReactionSystem`s
@@ -87,7 +87,9 @@ rx1 = Reaction(k,[B,C],[B,D], [2.5,1],[3.5, 2.5])
 rx2 = Reaction(2*k, [B], [D], [1], [2.5])
 rx3 = Reaction(2*k, [B], [D], [2.5], [2])
 @named mixedsys = ReactionSystem([rx1, rx2, rx3], t, [A, B, C, D], [k, b])
+mixedsys = complete(mixedsys)
 osys = convert(ODESystem, mixedsys; combinatoric_ratelaws = false)
+osys = complete(osys)
 ```
 Note, when using `convert(ODESystem, mixedsys; combinatoric_ratelaws=false)` the
 `combinatoric_ratelaws=false` parameter must be passed. This is also true when
@@ -127,6 +129,7 @@ t = default_t()
 rx1 = Reaction(β, [S, I], [I], [1,1], [2])
 rx2 = Reaction(ν, [I], [R])
 @named sir = ReactionSystem([rx1, rx2], t)
+sir = complete(sir)
 oprob = ODEProblem(sir, [], (0.0, 250.0))
 sol = solve(oprob, Tsit5())
 plot(sol)
@@ -162,7 +165,7 @@ Julia `Symbol`s corresponding to each variable/parameter to their values, or
 from ModelingToolkit symbolic variables/parameters to their values. Using
 `Symbol`s we have
 ```@example faq4
-using Catalyst, DifferentialEquations
+using Catalyst, OrdinaryDiffEq
 rn = @reaction_network begin
     α, S + I --> 2I
     β, I --> R
@@ -199,6 +202,7 @@ the second example, or one can use the `symmap_to_varmap` function to convert th
 `Symbol` mapping to a symbolic mapping. I.e. this works
 ```@example faq4
 osys = convert(ODESystem, rn)
+osys = complete(osys)
 
 # this works
 u0 = symmap_to_varmap(rn, [:S => 999.0, :I => 1.0, :R => 0.0])
@@ -221,6 +225,7 @@ rx1 = @reaction k, A --> 0
 rx2 = @reaction $f, 0 --> A
 eq = f ~ (1 + sin(t))
 @named rs = ReactionSystem([rx1, rx2, eq], t)
+rs = complete(rs)
 osys = convert(ODESystem, rs)
 ```
 In the final ODE model, `f` can be eliminated by using
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -89,7 +89,7 @@ Pkg.add("Catalyst")
 To solve Catalyst models and visualize solutions, it is also recommended to
 install DifferentialEquations.jl and Plots.jl
 ```julia
-Pkg.add("DifferentialEquations")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 ```
 
@@ -114,7 +114,7 @@ which in Jupyter notebooks will give the figure
 
 To generate and solve a mass action ODE version of the model we use
 ```@example ind1
-using DifferentialEquations
+using OrdinaryDiffEq
 p     = [:α => .1/1000, :β => .01]
 tspan = (0.0,250.0)
 u0    = [:S => 999.0, :I => 1.0, :R => 0.0]
diff --git a/docs/src/model_creation/compositional_modeling.md b/docs/src/model_creation/compositional_modeling.md
@@ -18,10 +18,10 @@ end
 Alternatively one can just build the `ReactionSystem` via the symbolic interface.
 ```@example ex0
 @parameters d
-@variable t
+t = default_t()
 @species X(t)
 rx = Reaction(d, [X], nothing)
-@named degradation_component = ReactionSystem([rs], t)
+@named degradation_component = ReactionSystem([rx], t)
 ```
 We can test whether a system is complete using the `ModelingToolkit.iscomplete` function:
 ```@example ex0
diff --git a/docs/src/model_creation/examples/hodgkin_huxley_equation.md b/docs/src/model_creation/examples/hodgkin_huxley_equation.md
@@ -13,7 +13,7 @@ cells such as neurons and muscle cells.
 We begin by importing some necessary packages.
 ```@example hh1
 using ModelingToolkit, Catalyst, NonlinearSolve
-using DifferentialEquations, Symbolics
+using OrdinaryDiffEq, Symbolics
 using Plots
 t = default_t()
 D = default_time_deriv()
@@ -77,6 +77,7 @@ I = I₀ * sin(2*pi*t/30)^2
 # get the gating variables to use in the equation for dV/dt
 @unpack m,n,h = hhrn
 
+Dₜ = default_time_deriv()
 eqs = [Dₜ(V) ~ -1/C * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + I/C]
 @named voltageode = ODESystem(eqs, t)
 nothing # hide
@@ -88,6 +89,7 @@ Finally, we add this ODE into the reaction model as
 
 ```@example hh1
 @named hhmodel = extend(voltageode, hhrn)
+hhmodel = complete(hhmodel)
 nothing # hide
 ```
 
diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md
@@ -7,7 +7,7 @@ use symbolic stoichiometries, and discuss several caveats to be aware of.
 Let's first consider a simple reversible reaction where the number of reactants
 is a parameter, and the number of products is the product of two parameters.
 ```@example s1
-using Catalyst, Latexify, DifferentialEquations, ModelingToolkit, Plots
+using Catalyst, Latexify, OrdinaryDiffEq, ModelingToolkit, Plots
 revsys = @reaction_network revsys begin
     k₊, m*A --> (m*n)*B
     k₋, B --> A
@@ -58,6 +58,7 @@ stoichiometries `(F,2*H,2)`.
 Let's now convert `revsys` to ODEs and look at the resulting equations:
 ```@example s1
 osys = convert(ODESystem, revsys)
+osys = complete(osys)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
@@ -88,6 +89,7 @@ converting to an `ODESystem`). For the previous example this gives the following
 (different) system of ODEs
 ```@example s1
 osys = convert(ODESystem, revsys; combinatoric_ratelaws = false)
+osys = complete(osys)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
@@ -139,7 +141,9 @@ The parameter `b` does not need to be explicitly declared in the
 
 We next convert our network to a jump process representation
 ```@example s1
+using JumpProcesses
 jsys = convert(JumpSystem, burstyrn; combinatoric_ratelaws = false)
+jsys = complete(jsys)
 equations(jsys)
 show(stdout, MIME"text/plain"(), equations(jsys)) # hide
 ```
diff --git a/docs/src/steady_state_functionality/bifurcation_diagrams.md b/docs/src/steady_state_functionality/bifurcation_diagrams.md
@@ -43,7 +43,7 @@ nothing # hide
 
 Finally, we compute our bifurcation diagram using:
 ```@example ex1
-bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
 nothing # hide
 ```
 Where `PALC()` designates that we wish to use the pseudo arclength continuation method to track our solution. The third argument (`2`) designates the maximum number of recursions when branches of branches are computed (branches appear as continuation encounters certain bifurcation points). For diagrams with highly branched structures (rare for many common small chemical reaction networks) this input is important. Finally, `bothside = true` designates that we wish to perform continuation on both sides of the initial point (which is typically the case). 
@@ -69,7 +69,7 @@ opt_newton = NewtonPar(tol = 1e-9, max_iterations = 1000)
 opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], 
                           dsmin = 0.001, dsmax = 0.01, max_steps = 1000,
                           newton_options = opt_newton)
-bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
 nothing # hide
 ```
 (however, in this case these additional settings have no significant effect on the result)
@@ -79,7 +79,7 @@ Let's consider the previous case, but instead compute the bifurcation diagram ov
 ```@example ex1
 p_span = (2.0, 15.0)
 opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
-bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
 plot(bif_dia; xguide = "k1", yguide = "X")
 ```
 Here, in the bistable region, we only see a single branch. The reason is that the continuation algorithm starts at our initial guess (here made at $k1 = 4.0$ for $(X,Y) = (5.0,2.0)$) and tracks the diagram from there. However, with the upper bound set at $k1=15.0$ the bifurcation diagram has a disjoint branch structure, preventing the full diagram from being computed by continuation alone. In this case it could be solved by increasing the bound from $k1=15.0$, however, this is not possible in all cases. In these cases, *deflation* can be used. This is described in the [BifurcationKit documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials2/#Snaking-computed-with-deflation).
@@ -103,7 +103,7 @@ bprob = BifurcationProblem(kinase_model, u_guess, p_start, :d; plot_var = :Xp, u
 
 p_span = (0.1, 10.0)
 opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
-bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
+bif_dia = bifurcationdiagram(bprob, PALC(), 2, opts_br; bothside = true)
 plot(bif_dia; xguide = "d", yguide = "Xp")
 ```
 This bifurcation diagram does not contain any interesting features (such as bifurcation points), and only shows how the steady state concentration of $Xp$ is reduced as $d$ increases. 
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
@@ -449,7 +449,7 @@ function remove_diffs(expr)
         return expr
     end
 end
-diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0.0 : expr)
+diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0 : expr)
 
 COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be converted to other system types. A ReactionSystem can be marked as complete using the `complete` function."
 
@@ -653,7 +653,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
 
     (remove_conserved !== nothing) &&
-        error("Catalyst does not support removing conserved species when converting to JumpSystems.")
+        throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems."))
 
     flatrs = Catalyst.flatten(rs)
     error_if_constraints(JumpSystem, flatrs)
diff --git a/test/dsl/dsl_basic_model_construction.jl b/test/dsl/dsl_basic_model_construction.jl
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
diff --git a/test/simulation_and_solving/simulate_ODEs.jl b/test/simulation_and_solving/simulate_ODEs.jl
diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl
diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl
diff --git a/test/test_networks.jl b/test/test_networks.jl
