You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Bifurcation diagrams describe how, for a dynamical system, the quantity and type of its steady states change as a parameter is varied[^1]. When using Catalyst-generated models, these can be computed with the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. Catalyst provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `ReactionSystem`s.
3
+
Bifurcation diagrams describe how, for a dynamical system, the quantity and type of its steady states change as a parameter is varied[^1]. When using Catalyst-generated models, these can be computed with the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. Catalyst provides a simple interface for creating BifurcationKit compatible [`BifurcationProblem`](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.BifurcationProblem)s from [`ReactionSystem`](@ref)s.
4
4
5
5
This tutorial briefly introduces how to use Catalyst with BifurcationKit through basic examples, with BifurcationKit.jl providing [a more extensive documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/). Especially for more complicated systems, where careful tuning of algorithm options might be required, reading the BifurcationKit documentation is recommended. Finally, BifurcationKit provides many additional features not described here, including [computation of periodic orbits](@ref bifurcationkit_periodic_orbits), [tracking of bifurcation points along secondary parameters](@ref bifurcationkit_codim2), and [bifurcation computations for PDEs](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials/#PDEs:-bifurcations-of-equilibria).
For this example, we will use a modified version of the model from Wilhelm (2009)[^2] (which
9
-
demonstrates a bistable switch as the parameter $k1$ is varied). We declare the model using Catalyst:
8
+
For this example, we will use a modified version of the [Wilhelm model](@ref basic_CRN_library_wilhelm) (which demonstrates a bistable switch as the parameter $k1$ is varied). We declare the model using Catalyst:
10
9
```@example ex1
11
10
using Catalyst
12
-
wilhelm_2009_model = @reaction_network begin
11
+
wilhelm_model = @reaction_network begin
13
12
k1, Y --> 2X
14
13
k2, 2X --> X + Y
15
14
k3, X + Y --> Y
@@ -21,23 +20,23 @@ Next we will create a `BifurcationProblem`. In addition to the `ReactionSystem`,
21
20
- The bifurcation parameter (the parameter which is varied in the bifurcation diagram).
22
21
- A full model parameter set. This includes the values of all non-bifurcation parameters, but also a value for the bifurcation parameter (which corresponds to the point in parameter space from which the computation of the bifurcation diagram starts).
23
22
- An initial guess of the steady state values of the system at the provided parameter set. Using this point as a starting guess for root finding, BifurcationKit calculates an initial steady state from which to compute the bifurcation diagram. Hence, this guess does not need to be very exact (but may be important if the system exhibits multistability for the initial parameter set).
24
-
- The species or statistic we wish to plot on the y-axis of the bifurcation diagram (a custom value can be provided by using the [`record_from_solution` argument](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#.-record*from*solution)).
23
+
- The species or statistic we wish to plot on the y-axis of the bifurcation diagram (customised properties can be recorded and plotted using the [`record_from_solution` argument](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#.-record*from*solution)).
25
24
26
25
We combine all this information to form a `BifurcationProblem`:
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:
36
+
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 two settings (`p_min` and `p_max`), which sets the minimum and maximum values over which the bifurcation parameter is varied, respectively. We wish to compute a bifurcation diagram over the interval $(2.0,20.0)$, and will use the following settings:
Where `PALC()` designates that we wish to use the pseudoarclength 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).
48
+
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).
50
49
51
50
We can plot our bifurcation diagram using the Plots.jl package:
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".
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:
60
59
-`p_min` and `p_max`: Set the interval over which the bifurcation diagram is computed (with the continuation stopping if it reaches these bounds).
61
60
-`dsmin` and `dsmax`: The minimum and maximum length of the continuation steps (in the bifurcation parameter's value).
61
+
-`ds`: The initial length of the continuation steps. This is especially important when `bothside = true`*is not* used, as teh sign of `ds` determines the direction from the initial point in which the continuation will proceed.
62
62
-`max_steps`: The maximum number of continuation steps. If a bifurcation diagram looks incomplete, try increasing this value.
63
63
-`newton_options`: Options for the Newton's method that BifurcationKit uses to find steady states. This can be created using `NewtonPar(tol = 1e-9, max_iterations = 100)` which here sets the tolerance (to `1e-9`) and the maximum number of newton iterations (to `100`).
64
64
65
65
The previous bifurcation diagram can be computed, with these various options specified, in the following way:
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).
84
+
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 = 8.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. Here we will describe a more general solution, the use of *deflation* to find multiple steady states.
86
85
86
+
We will not describe the details of how deflation works here (instead, please read [BifurcationKit's documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/DeflatedContinuation/)). In practise, in involves first creating a deflation continuation algorithm (`defalg`). This one, in turn, requires us to create [a deflation operator](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.DeflationOperator) and a perturbation function.
Next we compute our bifurcation diagram using `defalg` (instead of `PALC()` as previous) as our method. We also use a range of additional continuation parameter options to ensure a smooth tracking of the solution. Here we have
## [Systems with conservation laws](@id bifurcation_diagrams_cons_laws)
89
-
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.
106
+
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](@ref 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.
90
107
91
-
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$.
108
+
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$.
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.
110
127
111
128
Finally, for additional clarity, we reiterate the purpose of the two `u` arguments used:
112
-
-`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 fulfill the conserved concentrations provided in `u0`.
129
+
-`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`.
113
130
-`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.
114
131
115
132
116
133
---
117
-
## [Citation](@id bifurcation_kit_citation)
118
-
If you use this functionality in your research, please cite the following paper to support the author of the BifurcationKit package:
If you use BifurcationKit.jl for your work, we ask that you **cite** the following paper!! Open source development strongly depends on this. It is referenced on [HAL-Inria](https://hal.archives-ouvertes.fr/hal-02902346) with *bibtex* entry [CITATION.bib](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/CITATION.bib).
133
136
134
137
---
135
138
## References
136
-
[^1]: [Yuri A. Kuznetsov, *Elements of Applied Bifurcation Theory*, Springer (2023).](https://link.springer.com/book/10.1007/978-3-031-22007-4)
137
-
[^2]: [Thomas Wilhelm, *The smallest chemical reaction system with bistability*, BMC Systems Biology (2009).](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90)
139
+
[^1]: [Yuri A. Kuznetsov, *Elements of Applied Bifurcation Theory*, Springer (2023).](https://link.springer.com/book/10.1007/978-3-031-22007-4)
0 commit comments