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
Copy file name to clipboardExpand all lines: docs/src/steady_state_functionality/bifurcation_diagrams.md
+4-4Lines changed: 4 additions & 4 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -4,7 +4,7 @@ Bifurcation diagrams describe how, for a dynamical system, the quantity and type
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](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).
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
60
-`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
61
-`dsmin` and `dsmax`: The minimum and maximum length of the continuation steps (in the bifurcation parameter's value).
@@ -74,7 +74,7 @@ nothing # hide
74
74
```
75
75
(however, in this case these additional settings have no significant effect on the result)
76
76
77
-
## Bifurcation diagrams with disjoint branches
77
+
## [Bifurcation diagrams with disjoint branches](@id bifurcation_diagrams_disjoint_branches)
78
78
Let's consider the previous case, but instead compute the bifurcation diagram over the interval $(2.0,15.0)$:
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).
86
86
87
87
88
-
## Systems with conservation laws
88
+
## [Systems with conservation laws](@id bifurcation_diagrams_conservation_laws)
89
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.
90
90
91
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$.
We will first extract the equations corresponding to these from our model. Next, we will compute them using [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl). To generate our equations, we first convert our model to a `NonlinearSystem`. For each nullcline equation, we also have to replace the corresponding variable with a parameter (which can be varied to compute the full nullcline curve).
53
+
Here, each equation is an expression of two variables ($X$ and $Y$). We wish to plot $Y$ as a function of $X$. To do this, we will substitute the species variable $X$ with the parameter $Xpar$. This will enable us to carry out bifurcation analysis of the equations' solutions as $Xpar$ is varied
Next, we compute a `NonlinerSystem` corresponding to each nullcline.
59
+
To input these into BifurcationKit we need to convert these into `NonlinearSystem`s.
68
60
```@example nullcline_plotting
69
61
@named nc_X_sys = NonlinearSystem([nc_eq_X])
70
62
@named nc_Y_sys = NonlinearSystem([nc_eq_Y])
71
63
nc_X_sys = complete(nc_X_sys)
72
64
nc_Y_sys = complete(nc_Y_sys)
65
+
nothing # hide
73
66
```
74
-
Finally, for each of these, we can use BifurcationKit's `continuation` function to track the solution of the nullclines as the corresponding variable (the workflow is similar to when we used `bifurcationdiagram`[here](@ref bifurcation_diagrams)).
67
+
Finally, for out nullcline equations, we can use BifurcationKit's `continuation` function to track their solutions across all values of $Xpar$ (the workflow is similar to when we used `bifurcationdiagram`[here](@ref bifurcation_diagrams)).
75
68
```@example nullcline_plotting
76
69
using BifurcationKit
77
70
span = (0.0, 1.2)
78
-
function compute_nullcline(nc_sys, span, var_solve, pvar_bif)
Finally, we are ready to create our plot. We will plot the steady states using `scatter`. We use the `steady_state_stability` function to [compute steady state stabilities](@ref steady_state_stability) (and use this to determine how to plot the steady state markers).
81
+
We are ready to create our plot. We will plot the steady states using `scatter`. We use the `steady_state_stability` function to [compute steady state stabilities](@ref steady_state_stability) (and use this to determine how to plot the steady state markers).
89
82
```@example nullcline_plotting
90
83
using Plots
91
-
plot(nc_X; vars = (:x, :param), label = "dX/dt = 0", lw = 7, la = 0.7)
92
-
plot!(nc_Y; vars = (:param, :x), label = "dY/dt = 0", lw = 7, la = 0.7)
Here we can see how the steady states occur at the nullclines intersections.
106
99
100
+
!!! warn
101
+
BifurcationKit's `continuation` function will not detect disjoint branches, and above we can only use it because we know that each nullcline consists of a single branch. To handle disjoint nullclines, consider using [deflated continuation](@ref bifurcation_diagrams_disjoint_branches) (or possibly a package like [Contour.jl](https://github.com/JuliaGeometry/Contour.jl)).
102
+
107
103
## [Plotting system directions in phase space](@id nullcline_plotting_directions)
104
+
One useful property of nullclines is that the sign of $dX/dt$ will only switch whenever the solution crosses the $dX/dt=0$ nullcline. This mean that, within each region defined by the nullclines, the direction of the solution remains constant. Below we use this to, for each such region, plot arrows showing the solution's direction.
108
105
109
106
```@example nullcline_plotting
110
107
nlprob = NonlinearProblem(complete(nlsys), [X => 0.0, Y => 0.0], ps)
This also works as a form of simple stability analysis, where we can see how the solution moves *away* from the unstable steady state, and *to* the stable ones.
0 commit comments