|
| 1 | +# [Analysing model steady state properties with DynamicalSystems.jl](@id dynamical_systems) |
| 2 | +The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). For ODEs, methods primarily include those for analysing a system's steady state and stability properties. Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules). |
| 3 | + |
| 4 | +## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction) |
| 5 | +Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will (almost always) be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. |
| 6 | + |
| 7 | +DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create a `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` is unused, and its value is irrelevant. The value of `tspan` is not very important, however, it must provide long enough a time range (else basins will not be successfully found). |
| 8 | +```@example dynamical_systems_basins |
| 9 | +using Catalyst |
| 10 | +wilhelm_model = @reaction_network begin |
| 11 | + k1, Y --> 2X |
| 12 | + k2, 2X --> X + Y |
| 13 | + k3, X + Y --> Y |
| 14 | + k4, X --> 0 |
| 15 | +end |
| 16 | +
|
| 17 | +u0 = [:X => 0.0, :Y => .0] |
| 18 | +tspan = (0.0, 10.0) |
| 19 | +ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5] |
| 20 | +oprob = ODEProblem(wilhelm_model, u0, tspan, ps) |
| 21 | +nothing # hide |
| 22 | +``` |
| 23 | +Next, for any application of DynamicalSystems.jl, our `ODEProblem` must be converted into a so-called `CoupledODEs` structure. This is done by combining the ODE with the solver (and potential solver options) with which we wish to simulate it (just like when it is simulated using `solve`). Here, we will simply designate the `Tsit5` numeric solver (but provide no other options). |
| 24 | +```@example dynamical_systems_basins |
| 25 | +using DynamicalSystems, OrdinaryDiffEq |
| 26 | +ds = CoupledODEs(oprob, (alg = Tsit5(),)) |
| 27 | +nothing # hide |
| 28 | +``` |
| 29 | +We can now compute the basins of attraction. This is done by first designating a grid designate which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we first create a `AttractorsViaRecurrences` struct and then use that as input to the `basins_of_attraction` function. |
| 30 | +```@example dynamical_systems_basins |
| 31 | +# We provide one grid of values for each species. These are then bundled into a tuple. |
| 32 | +x_grid = 0.0:0.01:5.0 |
| 33 | +y_grid = 0.0:0.01:10.0 |
| 34 | +grid = (x_grid, y_grid) |
| 35 | +avr = AttractorsViaRecurrences(ds, grid) |
| 36 | +basins, attractors = basins_of_attraction(avr, grid) |
| 37 | +nothing # hide |
| 38 | +``` |
| 39 | +Here, `attractors` is a vector containing the attractors (in this case two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$). Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belong. |
| 40 | + |
| 41 | +DynamicalSystems.jl also provide a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), and alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a popular competitor to Plots.jl for normal plotting). |
| 42 | +```@example dynamical_systems_basins |
| 43 | +using CairoMakie |
| 44 | +heatmap_basins_attractors(grid, basins, attractors) |
| 45 | +``` |
| 46 | +Here, in addition to the basins of attraction, the system's three steady states are marked (the one at the intersection of the two basins is unstable). |
| 47 | + |
| 48 | +!!! warning |
| 49 | + Both Makie and Plots.jl exports a function called `plot`. Hence, if both these packages are imported into the same session, calls to `plot` must be prepended with the package one wishes to use (e.g. `Plots.plot(sol)`). |
| 50 | + |
| 51 | +More information on how to compute basins of attractions for ODEs using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/basins/). |
| 52 | + |
| 53 | +## [Computing Lyapunov exponents](@id dynamical_systems_lyapunov_exponents) |
| 54 | +[*Lyapunov exponents*](https://en.wikipedia.org/wiki/Lyapunov_exponent) are scalar values that can be computed for any state of an ODE. For an ODE with $n$ variables, for each state, a total of $n$ Lyapunov exponents can be computed (and these are collectively called the *Lyapunov spectrum*). Large Lyapunov exponents indicates that trajectories infinitesimally close to the given state diverge from each other. Conversely, small Lyapunov exponents suggests that trajectories converge to each others. |
| 55 | + |
| 56 | +While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviours in a point in phase space if that point have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not. |
| 57 | + |
| 58 | +First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviours. |
| 59 | +```@example dynamical_systems_lyapunov |
| 60 | +using Catalyst |
| 61 | +wr_model = @reaction_network begin |
| 62 | + k1, 2X --> 3X |
| 63 | + k2, X --> 2X |
| 64 | + k3, Z + 2X --> 2Z |
| 65 | + k4, Y + X --> 2Y |
| 66 | + k5, Y --> ∅ |
| 67 | + k6, 2Z --> ∅ |
| 68 | + k7, Z --> ∅ |
| 69 | +end |
| 70 | +``` |
| 71 | +We can simulate the model, noting that its behaviour seems chaotic. |
| 72 | +```@example dynamical_systems_lyapunov |
| 73 | +using OrdinaryDiffEq, Plots |
| 74 | +
|
| 75 | +u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5] |
| 76 | +tspan = (0.0, 100.0) |
| 77 | +p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7] |
| 78 | +
|
| 79 | +oprob = ODEProblem(wr_model, u0, tspan, p) |
| 80 | +sol = solve(oprob, Rodas5P()) |
| 81 | +plot(sol; idxs=(:X, :Y, :Z)) |
| 82 | +``` |
| 83 | +Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum (here, the `tspan` of the used `ODEProblem` have no impact on the computed Lyapunov spectrum). |
| 84 | +```@example dynamical_systems_lyapunov |
| 85 | +using DynamicalSystems |
| 86 | +ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff=false),)) |
| 87 | +nothing # hide |
| 88 | +``` |
| 89 | +Here, the `autodiff=false` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended). |
| 90 | +```@example dynamical_systems_lyapunov |
| 91 | +lyapunovspectrum(ds, 100) |
| 92 | +``` |
| 93 | +Here, the largest exponent is positive, suggesting that the model is chaotic in the point $(1.5, 1.5, 1.5)$. |
| 94 | + |
| 95 | +Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle: |
| 96 | +```@example dynamical_systems_lyapunov |
| 97 | +brusselator = @reaction_network begin |
| 98 | + A, ∅ --> X |
| 99 | + 1, 2X + Y --> 3X |
| 100 | + B, X --> Y |
| 101 | + 1, X --> ∅ |
| 102 | +end |
| 103 | +
|
| 104 | +u0_1 = [:X => 1.0, :Y => 1.0] |
| 105 | +u0_2 = [:X => 1.2, :Y => 1.0] |
| 106 | +tspan = (0., 25.) |
| 107 | +ps = [:A => 1.0, :B => 4.0] |
| 108 | +
|
| 109 | +oprob1 = ODEProblem(brusselator, u0_1, tspan, ps) |
| 110 | +oprob2 = ODEProblem(brusselator, u0_2, tspan, ps) |
| 111 | +osol1 = solve(oprob1, Tsit5()) |
| 112 | +osol2 = solve(oprob2, Tsit5()) |
| 113 | +plot(osol1; idxs = (:X, :Y)) |
| 114 | +plot!(osol2; idxs = (:X, :Y)) |
| 115 | +``` |
| 116 | +Next, we compute the Lyapunov spectrum at one of the initial conditions: |
| 117 | +```@example dynamical_systems_lyapunov |
| 118 | +ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff=false),)) |
| 119 | +lyapunovspectrum(ds, 100) |
| 120 | +``` |
| 121 | +Here, all Lyapunov exponents are negative, confirming that the brusselator (at this point) is non-chaotic. |
| 122 | + |
| 123 | +More details on how to compute Lyapunov exponents using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/lyapunovs/). A full overview of tools for analysing chaotic behaviours (using the "ChaosTools.jl subpackage) can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/). |
| 124 | + |
| 125 | + |
| 126 | +--- |
| 127 | +## [Citations](@id dynamical_systems_citations) |
| 128 | +If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following paper to support the author of the DynamicalSystems.jl package: |
| 129 | +``` |
| 130 | +@article{DynamicalSystems.jl-2018, |
| 131 | + doi = {10.21105/joss.00598}, |
| 132 | + url = {https://doi.org/10.21105/joss.00598}, |
| 133 | + year = {2018}, |
| 134 | + month = {mar}, |
| 135 | + volume = {3}, |
| 136 | + number = {23}, |
| 137 | + pages = {598}, |
| 138 | + author = {George Datseris}, |
| 139 | + title = {DynamicalSystems.jl: A Julia software library for chaos and nonlinear dynamics}, |
| 140 | + journal = {Journal of Open Source Software} |
| 141 | +} |
| 142 | +``` |
| 143 | + |
| 144 | + |
| 145 | +--- |
| 146 | +## References |
| 147 | +[^1]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7) |
| 148 | +[^2]: [S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Westview Press (1994).](http://users.uoa.gr/~pjioannou/nonlin/Strogatz,%20S.%20H.%20-%20Nonlinear%20Dynamics%20And%20Chaos.pdf) |
| 149 | +[^3]: [A. M. Lyapunov, *The general problem of the stability of motion*, International Journal of Control (1992).](https://www.tandfonline.com/doi/abs/10.1080/00207179208934253) |
0 commit comments