up

Author: TorkelE

docs/pages.jl

@@ -51,7 +51,8 @@ pages = Any[
         "steady_state_functionality/bifurcation_diagrams.md",
         "steady_state_functionality/dynamical_systems.md",
         "Examples" => Any[
-            "steady_state_functionality/examples/bifurcationkit_periodic_orbits.md"
+            "steady_state_functionality/examples/nullcline_plotting.md",
+            "steady_state_functionality/examples/bifurcationkit_periodic_orbits.md",
             "steady_state_functionality/examples/bifurcationkit_codim2.md"
         ]
     ],

docs/src/steady_state_functionality/examples/nullcline_plotting.md

@@ -1,5 +1,5 @@
 # [Plotting Nullclines and Steady States in Phase Space](@id nullcline_plotting)
-In this tutorial we will show how to extract a system's steady states and [nullclines](https://en.wikipedia.org/wiki/Nullcline), and how to plot these in [phase space](https://en.wikipedia.org/wiki/Phase_space). Generally, while nullclines are not directly needed for much analysis, plotting these can give some understanding of a system's steady state and stability properties. While nullclines can be "brute forced" (by solving their equations across grids of values), we will here use [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) to track their path exactly.
+In this tutorial we will show how to extract a system's steady states and [nullclines](https://en.wikipedia.org/wiki/Nullcline), and how to plot these in [phase space](https://en.wikipedia.org/wiki/Phase_space). Generally, while nullclines are not directly needed for much analysis, plotting these can give some understanding of a system's steady state and stability properties.
 
 For an ordinary differential equation
 ```math
@@ -13,20 +13,14 @@ For an ordinary differential equation
 the $i$'th nullcline is the surface along which $\frac{dx_i}{dt} = 0$, i.e. the implicit surface given by $f_i(x_1,\dots,x_n) = 0$. Nullclines are frequently used when visualizing the phase-planes of two-dimensional models (as these can be easily plotted).
 
 ## [Computing nullclines and steady states for a bistable switch](@id nullcline_plotting_computation)
-For our example we will use a simple bistable switch model, consisting of two species ($X$ and $Y$) which mutually inhibit each other through repressive Hill functions. We will create our model [programmatically](@ref programmatic_CRN_construction).
+For our example we will use a simple bistable switch model, consisting of two species ($X$ and $Y$) which mutually inhibit each other through repressive Hill functions. 
 ```@example nullcline_plotting
 using Catalyst
-t = default_t()
-@parameters v K n
-@species X(t) Y(t)
-rxs = [
-    Reaction(hillr(Y, v, K, n), [], [X]),
-    Reaction(hillr(X, v, K, n), [], [Y]),
-    Reaction(1.0, [X], []),
-    Reaction(1.0, [Y], []),
-]
-@named bs_switch = ReactionSystem(rxs, t)
-bs_switch = complete(bs_switch)
+bs_switch = @reaction_network begin
+    hillr(Y, v, K, n), 0 --> X
+    hillr(X, v, K, n), 0 --> Y
+    1.0, (X,Y) --> 0
+end
 ```
 
 Next, we compute the steady states [using homotopy continuation](@ref homotopy_continuation).
@@ -36,45 +30,23 @@ ps = [v => 1.0, K => 0.6, n => 4.0]
 sss = hc_steady_states(bs_switch, ps; show_progress = false)
 ```
 
-Finally, we will compute the nullclines. We will first extract our model's steady state equations (which we do by creating a `NonlinearSystem`).
-```@example nullcline_plotting
-nlsys = convert(NonlinearSystem, bs_switch)
-eqs = equations(nlsys)
-```
-Here, each equation is an expression of the 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
-```@example nullcline_plotting
-@parameters Xpar
-nc_eq_X = substitute(equations(nlsys)[2], Dict([X => Xpar]))
-nc_eq_Y = substitute(equations(nlsys)[1], Dict([X => Xpar]))
-```
-To input these into BifurcationKit we need to convert these into `NonlinearSystem`s.
+Finally, we will compute the nullclines. First we create a function which, for species values $(X,Y)$, returns the evaluation of the model's ODE's right-hand side.
 ```@example nullcline_plotting
-@named nc_X_sys = NonlinearSystem([nc_eq_X])
-@named nc_Y_sys = NonlinearSystem([nc_eq_Y])
-nc_X_sys = complete(nc_X_sys)
-nc_Y_sys = complete(nc_Y_sys)
-nothing # hide
-```
-Finally, for our 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)).
-```@example nullcline_plotting
-using BifurcationKit
-span = (0.0, 1.2)
-function compute_nullcline(nc_sys)
-    bprob = BifurcationProblem(nc_sys, [Y => 1.0], [ps; Xpar => 0.1], Xpar)
-    opts_br = ContinuationPar(p_min = span[1], p_max = span[2], dsmax = 0.01)
-    return continuation(bprob, PALC(), opts_br; bothside = true)
+nlprob = NonlinearProblem(complete(nlsys), [X => 0.0, Y => 0.0], ps)
+function get_XY(Xval, Yval)
+    prob = Catalyst.remake(nlprob; u0 = [X => Xval, Y => Yval])
+    return prob[equations(nlsys)[2].rhs], prob[equations(nlsys)[1].rhs]
 end
-nc_X = compute_nullcline(nc_X_sys)
-nc_Y = compute_nullcline(nc_Y_sys)
-nothing # hide
 ```
-
-We are now 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). For the nullclines, we plot their `.x` values (the variable we solve for, $Y$) against their `.param` values (the parameter which value we vary, $Xpar$).
+Next, we plot our nullclines by using Plot.jl's [`contour` function](https://docs.juliaplots.org/latest/series_types/contour/). Here, we plot the $0$ contour lines (which corresponds to the nullclines) for both the $X$ and $Y$ nullclines. 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).
 ```@example nullcline_plotting
 using Plots
-# Plot the nullclines.
-plot(nc_X.param, nc_X.x; label = "dX/dt = 0", lw = 7, la = 0.7)
-plot!(nc_Y.param, nc_Y.x; label = "dY/dt = 0", lw = 7, la = 0.7)
+# Plot the nullclines. Line labels added in separate `plot` commands (due to how the `contour` works).
+plt_mesh = 0:0.02:1.25
+contour(plt_mesh, plt_mesh, (x,y) -> get_XY(x,y)[1]; levels = [0.0], lw = 7, la = 0.7, color = 1, cbar=false)
+contour!(plt_mesh, plt_mesh, (x,y) -> get_XY(x,y)[2]; levels = [0.0], lw = 7, la = 0.7, color = 2, cbar=false)
+plot!([]; label = "dX/dt = 0", lw = 7, la = 0.7, color = 1)
+plot!([]; label = "dY/dt = 0", lw = 7, la = 0.7, color = 2)
 
 # Plot the steady states.
 for ss in sss
@@ -93,19 +65,12 @@ plot!(xlimit = span, ylimit = span, xguide = "X", yguide = "Y", legendfontsize =
 ```
 Here we can see how the steady states occur at the nullclines intersections.
 
-!!! warn
-    BifurcationKit's `continuation` function will not detect disjoint branches, and here we can only use it because we know that each nullcline consists of a single branch. To handle nullclines with disjoint branches, consider using [deflated continuation](@ref bifurcation_diagrams_disjoint_branches) (or possibly a package like [Contour.jl](https://github.com/JuliaGeometry/Contour.jl)).
+!!! note
+    Here we use an inherent Plots function to plot the nullclines. However, there are also specialised packages for these kinds of plots, such as [ImplicitPlots.jl](https://github.com/saschatimme/ImplicitPlots.jl).
 
 ## [Plotting system directions in phase space](@id nullcline_plotting_directions)
 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 means 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.
 ```@example nullcline_plotting
-# Create a function which evaluates the ODE at a point in phase space.
-nlprob = NonlinearProblem(complete(nlsys), [X => 0.0, Y => 0.0], ps)
-function get_XY(Xval, Yval)
-    prob = Catalyst.remake(nlprob; u0 = [X => Xval, Y => Yval])
-    return prob[equations(nlsys)[2].rhs], prob[equations(nlsys)[1].rhs]
-end
-
 # Creates a function for plotting the ODE's direction at a point in phase space.
 function plot_xy_arrow!(Xval, Yval)
     dX, dY = get_XY(Xval, Yval)