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/examples/bifurcationkit_codim2.md
+16-16Lines changed: 16 additions & 16 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -1,9 +1,9 @@
1
-
# # [Tracking Bifurcation Point w.r.t. Secondary Parameters using BifurcationKit](@idbifurcationkit_bifpoint_continuation)
1
+
# [Tracking Bifurcation Point w.r.t. Secondary Parameters using BifurcationKit.jl](@idbifurcationkit_codim2)
2
2
Previously, we have shown how to [compute bifurcation diagrams](@ref bifurcation_diagrams) using [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl). In this example, we will show how, after computing the initial diagram, we can track how the position of a bifurcation point moves as a secondary parameter is changed (so-called codimensional 2 bifurcation analysis). More information on how to track bifurcation points along secondary parameters can be found in the [BifurcationKit documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/tutorials/ode/tutorialCO/#CO-oxidation-(codim-2)).
3
3
4
-
## [Computing the bifurcation diagram for the Repressilator](@idbifurcationkit_bifpoint_continuation_bifdia)
4
+
## [Computing the bifurcation diagram for the Repressilator](@idbifurcationkit_codim2_bifdia)
5
5
We will first compute the bifurcation diagram, using the same approach as in the [corresponding tutorial](@ref bifurcation_diagrams). For this example, we will use the oscillating [Repressilator](@ref basic_CRN_library_repressilator) model.
6
-
```@examplebifurcationkit_bifpoint_continuation
6
+
```@examplebifurcationkit_codim2
7
7
using Catalyst
8
8
repressilator = @reaction_network begin
9
9
hillr(Z,v,K,n), ∅ --> X
@@ -13,7 +13,7 @@ repressilator = @reaction_network begin
13
13
end
14
14
```
15
15
Next, we create a `BifurcationProblem` for our model. We will compute the bifurcation diagram with respect to the parameter $v$, and plot the species $X$ in the diagram.
We note that for small values of $v$ the system's single steady state is stable (where the line is thicker). After a [Hopf](https://en.wikipedia.org/wiki/Hopf_bifurcation) bifurcation (the red point), the state turns unstable (where the line is thinner). For chemical reaction networks (which mostly are well-behaved) a single unstable steady state typically corresponds to an oscillation. We can confirm that the system oscillates in the unstable region (while it reaches a stable steady state in the stable region) using simulations:
plot(plot(sol_nosc; title = "v = 5"), plot(sol_osc; title = "v = 15"), size = (1000,400), lw = 4)
48
48
```
49
49
50
-
## [Tracking the bifurcation point w.r.t. a second parameter](@idbifurcationkit_bifpoint_continuation_codim2_2ndpar_cont)
50
+
## [Tracking the bifurcation point w.r.t. a second parameter](@idbifurcationkit_codim2_2ndpar_cont)
51
51
Next, we will investigate how the Hopf bifurcation point moves (in $v$-$X$ space) as a second parameter ($K$) is changed. To do this we will use BifurcationKit.jl's [`continuation` function](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/library/#BifurcationKit.continuation) (the [`bifurcationdiagram` function](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/library/#BifurcationKit.bifurcationdiagram), which we previously have used, works by calling `continuation` recursively). We will call it on the Hopf bifurcation point. First, we will designate the parameter we wish to change, as well as which interval to change it. For this, we compute the parameter index (`K_idx`) that BifurcationKit will see (a smoother interface for this will hopefully be added in the future). We also create a new `ContinuationPar`, to which we add our $K$ interval.
Now we can compute the continuation of the Hopf bifurcation. First we must extract a branch from our bifurcation diagram (using [`get_branch`](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.get_branch)), as `continuation` cannot work on bifurcation diagrams directly (our bifurcation diagram consists of a single branch, which we here extract). Next, we can call `continuation` on it, designating that we wish to perform continuation from the second point on the branch (which corresponds to the Hopf bifurcation, the first point is one of the branch's two endpoints).
We can now plot how the position of the bifurcation point changes with $K$. Here, we use `vars = (:p1, :x)` to designate that we wish to plot (across the continuation branch) the plotting variable ($X$, which we designated when we created our `BifurcationProblem`) against the first parameter ($v$).
65
-
```@examplebifurcationkit_bifpoint_continuation
65
+
```@examplebifurcationkit_codim2
66
66
plot(bifdia; branchlabel = "Continuation of steady state w.r.t. v, (K = 15)", linewidthstable = 6,
@@ -71,7 +71,7 @@ plot!(cont_hopf; vars = (:p1, :x), xlimit = v_span, xguide = bif_par, yguide = p
71
71
In this case we cannot see directly which part of the $K$ continuation branch corresponds to low values, however, for low $K$ the Hopf bifurcation occurs for much lower values of $v$ (and corresponds to lower steady state values of $X$). We can check this by e.g. re-computing the Hopf branch for `K_span = (0.01, 20.0)` and see that the rightmost part of the branch is shortened.
72
72
73
73
We can confirm that the new line corresponds to the Hopf Bifurcation point by recomputing the initial bifurcation diagram, but for a lower $K$ value.
Here we see that the Hopf bifurcation point of the new diagram also lies on the Hopf continuation line.
82
82
83
83
Finally, we have already noted that the Hopf bifurcation splits parameter space into one part where the system oscillates and one where it doesn't. Previously we plotted the Hopf continuation in $v$-$X$ space, however, it is also possible to plot it in $v$-$K$ space using the `vars = (:p1, :p2)` argument:
Next, we colour parameter space according to whether the steady state is stable (blue) or unstable (red). We also mark two sample values (one in each region).
91
-
```@examplebifurcationkit_bifpoint_continuation
91
+
```@examplebifurcationkit_codim2
92
92
sample1 = (15.0, 10.0)
93
93
sample2 = (5.0, 15.0)
94
94
plot(cont_hopf; vars = (:p1, :p2), fillrange = ylimit[2])
Finally, we can perform one simulation using each of the parameter samples, confirming that one corresponds to an oscillation, while the other one does not.
0 commit comments