up

Project.toml

@@ -32,11 +32,13 @@ CatalystBifurcationKitExtension = "BifurcationKit"
 CatalystHomotopyContinuationExtension = "HomotopyContinuation"
 
 [compat]
+BifurcationKit = "0.3"
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
 Graphs = "1.4"
 JumpProcesses = "9.3.2"
+HomotopyContinuation = "2.9"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"

docs/src/catalyst_applications/bifurcation_diagrams.md

@@ -6,7 +6,7 @@ This tutorial briefly introduces how to use Catalyst with BifurcationKit through
 
 ## Basic example
 For this example, we will use a modified version of the model from Wilhelm (2009)[^2] (which
-demonstrates a bistable switch as the parameter *k1* is varied). We declare the model using Catalyst:
+demonstrates a bistable switch as the parameter $k1$ is varied). We declare the model using Catalyst:
 ```@example ex1
 using Catalyst
 wilhelm_2009_model = @reaction_network begin
@@ -34,7 +34,7 @@ bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_v
 nothing # hide
 ```
 
-BifurcationKit computes bifurcation diagrams using the `bifurcationdiagram` function. From an initial point in the diagram, it tracks the solution (using a continuation algorithm) until the entire diagram is computed (BifurcationKit's continuation can be used for other purposes, however, this tutorial focuses on bifurcation diagram computation). The continuation settings are provided in a `ContinuationPar` structure. In this example, we will only specify three settings, `p_min` and `p_max` (which sets the minimum and maximum values over which the bifurcation parameter is varied) and `max_steps` (the maximum number of continuation steps to take as the bifurcation diagram is tracked). We wish to compute a bifurcation diagram over the interval *(2.0,20.0)*, and will use the following settings:
+BifurcationKit computes bifurcation diagrams using the `bifurcationdiagram` function. From an initial point in the diagram, it tracks the solution (using a continuation algorithm) until the entire diagram is computed (BifurcationKit's continuation can be used for other purposes, however, this tutorial focuses on bifurcation diagram computation). The continuation settings are provided in a `ContinuationPar` structure. In this example, we will only specify three settings, `p_min` and `p_max` (which sets the minimum and maximum values over which the bifurcation parameter is varied) and `max_steps` (the maximum number of continuation steps to take as the bifurcation diagram is tracked). We wish to compute a bifurcation diagram over the interval $(2.0,20.0)$, and will use the following settings:
 ```@example ex1
 p_span = (2.0, 20.0)
 opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps=1000)
@@ -53,7 +53,7 @@ We can plot our bifurcation diagram using the Plots.jl package:
 using Plots
 plot(bif_dia; xguide="k1", yguide="X")
 ```
-Here, the steady state concentration of *X* is shown as a function of *k1*'s value. Stable steady states are shown with thick lines, unstable ones with thin lines. The two fold bifurcation points are marked with "bp".
+Here, the steady state concentration of $X$ is shown as a function of $k1$'s value. Stable steady states are shown with thick lines, unstable ones with thin lines. The two [fold bifurcation points](https://en.wikipedia.org/wiki/Saddle-node_bifurcation) are marked with "bp".
 
 ## Additional `ContinuationPar` options
 Most of the options required by the `bifurcationdiagram` function are provided through the `ContinuationPar` structure. For full details, please read the [BifurcationKit documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/library/#BifurcationKit.ContinuationPar). However, a few common options, and how they affect the continuation computation, are described here:
@@ -75,18 +75,18 @@ nothing # hide
 (however, in this case these additional settings have no significant effect on the result)
 
 ## Bifurcation diagrams with disjoint branches
-Let's consider the previous case, but instead compute the bifurcation diagram over the interval *(2.0, 15.0)*:
+Let's consider the previous case, but instead compute the bifurcation diagram over the interval $(2.0,15.0)$:
 ```@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)
 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).
+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).
 
 
 ## Systems with conservation laws
-Some systems are under-determined determined at steady-state, so that for a given parameter set they have an infinite number of possible steady state solutions, preventing bifurcation diagrams from being computed. Similar to when we [compute steady states for fixed parameter values](@ref homotopy_continuation_conservation_laws), we can utilise Catalyst's ability to detect and eliminate conservation laws to resolve this issue. This requires us to provide information of the species concentrations at which we wish to compute the bifurcation diagram (to determine the values of conserved quantities). These are provided to the `BifurcationProblem` using the `u0` argument.
+Some systems are under-determined at steady-state, so that for a given parameter set they have an infinite number of possible steady state solutions, preventing bifurcation diagrams from being computed. Similar to when we [compute steady states for fixed parameter values](@ref homotopy_continuation_conservation_laws), we can utilise Catalyst's ability to detect and eliminate conservation laws to resolve this issue. This requires us to provide information of the species concentrations at which we wish to compute the bifurcation diagram (to determine the values of conserved quantities). These are provided to the `BifurcationProblem` using the `u0` argument.
 
 To illustrate this, we will create a simple model of a kinase that is produced and degraded (at rates $p$ and $d$). The kinase facilitates the phosphorylation of a protein ($X$), which is dephosphorylated at a constant rate. For this system, we will compute a bifurcation diagram, showing how the concentration of the phosphorylated protein ($Xp$) depends on the degradation rate of the kinase ($d$). We will set the total amount of protein ($X+Xp$) to $1.0$.
 ```@example ex2
@@ -106,9 +106,11 @@ 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)
 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. For this example, we will note two facts:
-- When providing the concentrations for computing the conserved quantities (in `u0`), we only have to designate the concentrations of species that are actually involved in conservation laws. For larger systems, determining which one are may, however, be difficult. In this case, it might still be wise to provide concentrations for all species.
-- The steady state guess in `u_guess` does not actually have to fulfil the conserved concentrations provided in `u0`.
+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. 
+
+Finally, for additional clarity, we reiterate the purpose of the two `u` arguments used:
+- `u_guess`: A guess of the initial steady states (which BifurcationKit uses to find its starting point). Typically, most trivial guesses work (e.g. setting all species concentrations to `1.0`). `u_guess` *does not* have to fulfil the conserved concentrations provided in `u0`.
+- `u0`: Used to compute the concentrations of any conserved quantities (e.g. in our example $X + Xp = 1.0$). Technically, values are only required for species that are involved in conservation laws (in our case we do not need to provide a value for $K$). However, sometimes determining which species are actually involved in conservation laws can be difficult, and it might be easier to simply provide concentrations for all species.
 
 !!! note
     Computing bifurcation diagrams for [hierarchical models created by composing multiple systems](@ref compositional_modeling), that also contain conservations laws, is currently not supported. For these, please [flatten](@ref api_network_extension_and_modification) the system first.

src/miscellaneous.jl

@@ -1,13 +1 @@
 # File containing code which I am currently unsure which file in which it belongs. Long-term functions here should be moved elsewhere.
-
-
-# If u0s are not given while conservation laws are present, throws an error.
-# If the system is nested, the function is skipped (conservation laws cannot be computed for nested systems).
-# Used in HomotopyContinuation and BifurcationKit extensions.
-function conservationlaw_errorcheck(rs, pre_varmap)
-    isempty(get_systems(rs)) || return
-    vars_with_vals = Set(p[1] for p in pre_varmap)
-    any(s -> s in vars_with_vals, species(rs)) && return
-    isempty(conservedequations(rs)) || 
-        error("The system has conservation laws but initial conditions were not provided for some species.")
-end

src/networkapi.jl

@@ -1277,6 +1277,18 @@ Compute conserved quantities for a system with the given conservation laws.
 """
 conservedquantities(state, cons_laws) = cons_laws * state
 
+
+# If u0s are not given while conservation laws are present, throws an error.
+# If the system is nested, the function is skipped (conservation laws cannot be computed for nested systems).
+# Used in HomotopyContinuation and BifurcationKit extensions.
+function conservationlaw_errorcheck(rs, pre_varmap)
+    isempty(get_systems(rs)) || return
+    vars_with_vals = Set(p[1] for p in pre_varmap)
+    any(s -> s in vars_with_vals, species(rs)) && return
+    isempty(conservedequations(rs)) || 
+        error("The system has conservation laws but initial conditions were not provided for some species.")
+end
+
 ######################## reaction network operators #######################
 
 """

test/runtests.jl

@@ -21,9 +21,7 @@ using SafeTestsets
     # Runs various miscellaneous tests.
     @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end
     @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end
-    if VERSION >= v"1.9"
-        @time @safetestset "NonlinearProblems and Steady State Solving" begin include("miscellaneous_tests/nonlinear_solve.jl") end
-    end
+    @time @safetestset "NonlinearProblems and Steady State Solving" begin include("miscellaneous_tests/nonlinear_solve.jl") end
     @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end
     @time @safetestset "Compound species" begin include("miscellaneous_tests/compound_macro.jl") end
     @time @safetestset "Reaction balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end
@@ -54,8 +52,7 @@ using SafeTestsets
     end
 
     ### Tests extensions. ###
-    if VERSION >= v"1.9"
-        @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
-        @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
-    end
+    @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
+    @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
+
 end # @time