Skip to content

Commit 419edfc

Browse files
TorkelETorkelE
committed
add documentation
1 parent 1202c04 commit 419edfc

File tree

2 files changed

+89
-88
lines changed

2 files changed

+89
-88
lines changed

‎docs/src/catalyst_applications/bifurcation_diagrams.md

Lines changed: 88 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -1,115 +1,114 @@
11
# [Bifurcation Diagrams](@id bifurcation_diagrams)
2-
Bifurcation diagrams can be produced for models generated by Catalyst through the
3-
use of the [BifurcationKit.jl](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/)
4-
package. This tutorial gives a simple example of how to create such a
5-
bifurcation diagram.
62

7-
!!! note
8-
Catalyst 13.0 and up require at least BifurcationKit v0.2.4.
3+
Bifurcation diagrams describes how, for a dynamic system, the quantity and quality of its steady states changes with a parameter's value[^1]. These can be computed through the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. Catalyst provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `ReactionSystem`s. If you use this in your research, please [cite the BifurcationKit.jl](@ref bifurcation_kit_citation) and [Catalyst.jl](@ref catalyst_citation) publications.
94

10-
First, we declare our reaction model. For this example we will use a bistable
11-
switch, but one which also contains a Hopf bifurcation.
5+
This tutorial briefly introduces how to use Catalyst in conjecture 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](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/), [tracking of bifurcation points along secondary parameters](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/branchswitching/), and [bifurcation computations for PDEs](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials/#PDEs:-bifurcations-of-equilibria).
6+
7+
## Basic example
8+
9+
For this example, we will use a modified version of the model from Wilhelm (2009)[^2] (which
10+
demonstrates a bistable switch as the parameter *k1* is varied). We declare the model using Catalyst:
1211
```@example ex1
1312
using Catalyst
14-
rn = @reaction_network begin
15-
(v0 + v*(S * X)^n / ((S*X)^n + (D*A)^n + K^n), d), ∅ ↔ X
16-
(X/τ, 1/τ), ∅ ↔ A
13+
wilhelm_2009_model = @reaction_network begin
14+
k1, Y --> 2X
15+
k2, 2X --> X + Y
16+
k3, X + Y --> Y
17+
k4, X --> 0
18+
k5, 0 --> X
1719
end
1820
```
19-
Next, we specify the system parameters for which we wish to plot the bifurcation
20-
diagram. We also set the parameter we wish to vary in our bifurcation diagram,
21-
as well as the interval to vary it over. Finally, we set which variable we wish
22-
to plot the steady state values of in the bifurcation plot.
21+
Next we will create a `BifurcationProblem`. In addition to the `REactionSystem`, we need to provide:
22+
- The bifurcation parameter (the parameter which is varied in the bifurcation diagram).
23+
- A full model parameter set. This includes the values of all non-bifurcation parameters, but also a value for the bifurcation parameter (which correspond to the point in parameter space from which the computation of the bifurcation diagram starts).
24+
- An initial guess of the steady state values of the system at the provided parameter set. Using this point as a start, BifurcationKit uses Newton's method to find 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).
25+
- The species which concentration we wish to plot on the y-axis of the bifurcation diagram (alternatively, a custom value can be provided by using the [`record_from_solution` argument](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#.-record*from*solution)).
26+
27+
We combines all this information to a `BifurcationProblem`:
2328
```@example ex1
24-
p = Dict(:S => 1., :D => 9., :τ => 1000., :v0 => 0.01,
25-
:v => 2., :K => 20., :n => 3, :d => 0.05)
26-
bif_par = :S # bifurcation parameter
27-
p_span = (0.1, 20.) # interval to vary S over
28-
plot_var = :X # we will plot X vs S
29-
nothing # hide
29+
using BifurcationKit
30+
bif_par = :k1
31+
u_guess = [:X => 5.0, :Y => 2.0]
32+
p_start = [:k1 => 4.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.5, :k5 => 1.25]
33+
plot_var = :X
34+
bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var=plot_var)
35+
nothing # hide
3036
```
31-
When creating a bifurcation diagram, we typically start at some point in
32-
parameter-space. We will simply select the beginning of the interval over which
33-
we wish to compute the bifurcation diagram, `p_span[1]`. We thus create a
34-
modified parameter set where `S = 0.1`. For this parameter set, we also guess
35-
the steady state of the system. While a good estimate could be provided through
36-
an ODE simulation, BifurcationKit does not require the guess to be very
37-
accurate.
37+
38+
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). If we wish to compute a bifurcation diagram over the interval *(2.0,20.0)* we use the following settings:
3839
```@example ex1
39-
p_bstart = copy(p)
40-
p_bstart[bif_par] = p_span[1]
41-
u0 = [:X => 1.0, :A => 1.0]
42-
nothing # hide
40+
p_span = (2.0, 20.0)
41+
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps=1000)
42+
nothing # hide
4343
```
44-
Finally, we extract the ODE derivative function and its jacobian in a form that
45-
BifurcationKit can use:
44+
45+
Finally, we compute our bifurcation diagram using:
4646
```@example ex1
47-
oprob = ODEProblem(rn, u0, (0.0, 0.0), p_bstart; jac = true)
48-
F = (u,p) -> oprob.f(u, p, 0)
49-
J = (u,p) -> oprob.f.jac(u, p, 0)
50-
nothing # hide
47+
bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
48+
nothing # hide
5149
```
50+
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 encounter certain bifurcation points). For diagrams with highly branches structure (rare for most standard 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).
5251

53-
In creating an `ODEProblem` an ordering is chosen for the initial condition and
54-
parameters, and regular `Float64` vectors of their numerical values are created
55-
as `oprob.u0` and `oprob.p` respectively. BifurcationKit needs to know the index
56-
in `oprob.p` of our bifurcation parameter, `:S`, and the index in `oprob.u0` of
57-
the variable we wish to plot, `:X`. We calculate these as
52+
We can plot our bifurcation diagram using the Plots.jl package:
5853
```@example ex1
59-
# get S and X as symbolic variables
60-
@unpack S, X = rn
61-
62-
# find their indices in oprob.p and oprob.u0 respectively
63-
bif_idx = findfirst(isequal(S), parameters(rn))
64-
plot_idx = findfirst(isequal(X), species(rn))
65-
nothing # hide
54+
using Plots
55+
plot(bif_dia; xguide="k1", yguide="X")
6656
```
57+
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".
6758

68-
Now, we load the required packages to create and plot the bifurcation diagram.
69-
We also bundle the information we have compiled so far into a
70-
`BifurcationProblem`.
71-
```@example ex1
72-
using BifurcationKit, Plots, LinearAlgebra, Setfield
59+
## Additional `ContinuationPar` options
60+
Most of the options required by the `bifurcationdiagram` function is not provided directly, but instead 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:
61+
- `p_min` and `p_max`: Sets the interval over which the bifurcation diagram is computed (with the continuation stopping if it reaches these bounds).
62+
- `dsmin` and `dsmax`: The minimum and maximum length of the continuation steps (in the bifurcation parameter's value).
63+
- `max_steps`: The maximum number of continuation steps. If a bifurcation diagram looks incomplete, try increasing this value.
64+
- `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`).
7365

74-
bprob = BifurcationProblem(F, oprob.u0, oprob.p, (@lens _[bif_idx]);
75-
record_from_solution = (x, p) -> x[plot_idx], J = J)
76-
nothing # hide
77-
```
78-
Next, we need to specify the input options for the pseudo-arclength continuation
79-
method (PACM) which produces the diagram.
66+
The previous bifurcation diagram can be computed, with these various options specified, in the following way:
8067
```@example ex1
81-
bopts = ContinuationPar(dsmax = 0.05, # Max arclength in PACM.
82-
dsmin = 1e-4, # Min arclength in PACM.
83-
ds = 0.001, # Initial (positive) arclength in PACM.
84-
max_steps = 100000, # Max number of steps.
85-
p_min = p_span[1], # Min p-val (if hit, the method stops).
86-
p_max = p_span[2], # Max p-val (if hit, the method stops).
87-
detect_bifurcation = 3) # Value in {0,1,2,3}
88-
nothing # hide
68+
p_span = (2.0, 20.0)
69+
opt_newton = NewtonPar(tol = 1e-9, max_iterations = 1000)
70+
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2],
71+
dsmin = 0.001, dsmax = 0.01, max_steps = 1000,
72+
newton_options = opt_newton)
73+
bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
74+
nothing # hide
8975
```
90-
Here `detectBifurcation` determines to what extent bifurcation points are
91-
detected and how accurately their values are determined. Three indicates to use the most
92-
accurate method for calculating their values.
9376

94-
We are now ready to compute the bifurcation diagram:
77+
## Bifurcation diagrams with disjoint branches
78+
Let's consider the previous case, but instead compute the bifurcation diagram over the interval *(2.0, 15.0)*:
9579
```@example ex1
96-
bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> bopts)
97-
nothing # hide
98-
```
99-
Finally, we can plot it:
100-
```@example ex1
101-
plot(bf; xlabel = string(bif_par), ylabel = string(plot_var))
80+
p_span = (2.0, 15.0)
81+
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
82+
bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
83+
plot(bif_dia; xguide="k1", yguide="X")
10284
```
85+
Here, in the bistable region, we only see a single branch. The reason is that the continuation algorithm start at our initial guess (here made at *k1 = 4.0* for *(X,Y) = (5.0,2.0)*) and tracks diagram from there. However, with the upper bound set at *k1=15.0* the bifurcation diagram have 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 of *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).
10386

104-
Here, the Hopf bifurcation is marked with a red dot and the fold bifurcations
105-
with blue dots. The region with a thinner line width corresponds to unstable
106-
steady states.
10787

108-
This tutorial demonstrated how to make a simple bifurcation diagram where all
109-
branches are connected. However, BifurcationKit.jl is a very powerful package
110-
capable of a lot more. For more details, please see that package's
111-
documentation:
112-
[BifurcationKit.jl](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/).
88+
## Systems with conservation laws
89+
Some systems are under-determined, and for a given parameter set they have an infinite number of possible steady states, preventing bifurcation diagrams from being computed. Similar to when we [compute single steady states](@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. These are provided to the `BifurcationProblem` using the `u0` argument.
90+
91+
The 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 phosphoryalted protein (*Xp*) depends on teh degradation rate of the kinase (*d*). We will set the total amount of protein (*X+Xp*) to *1.0*.
92+
```@example ex2
93+
using BifurcationKit, Catalyst, Plots
94+
kinase_model = @reaction_network begin
95+
(p, d), 0 <--> K
96+
(K*kP,kD), X <--> Xp
97+
end
98+
99+
u_guess = [:K => 1.0, :X => 1.0, :Xp => 1.0]
100+
p_start = [:p => 1.0, :d => 0.5, :kP => 2.0, :kD => 5.0]
101+
u0 = [:X => 1.0, :Xp => 0.0]
102+
bprob = BifurcationProblem(kinase_model, u_guess, p_start, :d; plot_var=:Xp, u0=u0)
103+
104+
p_span = (0.1, 10.0)
105+
opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
106+
bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
107+
plot(bif_dia; xguide="d", yguide="Xp")
108+
```
109+
This bifurcation diagram does not contain any interesting features (such as bifurcation points), but only shows how the steady state concentration of *Xp* is reduced as *d* increases. For this example, we will note two additional facts:
110+
- 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 be wise to provide concentrations for all species.
111+
- The steady state guess in `u_guess` does not actually have to fulfil the conserved concentrations provided in `u0`.
113112

114113
---
115114
## [Citation](@id bifurcation_kit_citation)
@@ -131,4 +130,5 @@ If you use this functionality in your research, please cite the following paper
131130

132131
---
133132
## References
134-
[^1]: [Yuri A. Kuznetsov, *Elements of Applied Bifurcation Theory*, Springer (2023).](https://link.springer.com/book/10.1007/978-3-031-22007-4)
133+
[^1]: [Yuri A. Kuznetsov, *Elements of Applied Bifurcation Theory*, Springer (2023).](https://link.springer.com/book/10.1007/978-3-031-22007-4)
134+
[^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)

‎test/extensions/bifurcation_kit.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ end
5252
# Bistable switch.
5353
# Checks that the same bifurcation problem is created as for BifurcationKit.
5454
# Checks with Symbolics as bifurcation and plot vars.
55+
# Tries setting `jac=false`.
5556
let
5657
# Creates BifurcationProblem via Catalyst.
5758
bistable_switch = @reaction_network begin

0 commit comments

Comments
 (0)