init

Author: TorkelE

docs/pages.jl

@@ -42,7 +42,10 @@ pages = Any[
         "steady_state_functionality/nonlinear_solve.md",
         # "steady_state_functionality/steady_state_stability_computation.md",
         "steady_state_functionality/bifurcation_diagrams.md",
-        "steady_state_functionality/dynamical_systems.md"
+        "steady_state_functionality/dynamical_systems.md",
+        "Examples" => Any[
+            "steady_state_functionality/examples/nullcline_plotting.md"
+        ]
     ],
     "Inverse problems" => Any[
         "inverse_problems/petab_ode_param_fitting.md",

docs/src/steady_state_functionality/examples/nullcline_plotting.md

@@ -0,0 +1,91 @@
+# [Plotting Nullclines and Steady States in Phase Space](@id nullcline_plotting)
+
+Introduction text.
+
+For our example we will use a simple bistable switch model, consisting of two species ($X$ and $Y$) which mutually inhibits each other through repressive hill functions. We will create our model [programmatically](@ref programmatic_CRN_construction).
+```@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)
+```
+Next, we compute the steady states [using homotopy continuation](@ref homotopy_continuation).
+```@example nullcline_plotting
+import HomotopyContinuation
+ps = [v => 1.0, K => 0.6, n => 4.0]
+sss = hc_steady_states(bs_switch, ps; show_progress = false)
+```
+Next we wish to compute the nullclines. 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).
+```@example nullcline_plotting
+@parameters Xval Yval
+nlsys = convert(NonlinearSystem, bs_switch)
+nc_eq_X = substitute(equations(nlsys)[1], Dict([X => Xval]))
+nc_eq_Y = substitute(equations(nlsys)[2], Dict([Y => Yval]))
+```
+Next, we compute a `NonlinerSystem` corresponding to each nullcline.
+```@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)
+```
+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)).
+```@example nullcline_plotting
+using BifurcationKit
+span = (0.0, 1.2)
+function compute_nullcline(nc_sys, span, var_solve, pvar_bif)
+    bprob = BifurcationProblem(nc_sys, [var_solve => 1.0], [ps; pvar_bif => 1.0], pvar_bif; plot_var = var_solve)
+    opts_br = ContinuationPar(p_min = span[1], p_max = span[2], dsmax = 0.01) # `dsmax = 0.01` ensures a smooth plot.
+    return continuation(bprob, PALC(), opts_br; bothside = true)
+end
+nc_X = compute_nullcline(nc_X_sys, span, Y, Xval)
+nc_Y = compute_nullcline(nc_Y_sys, span, X, Yval)
+nothing # hide
+```
+
+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).
+```@example nullcline_plotting
+using Plots
+plot(nc_X; vars = (:x, :param), label = "dX/dt = 0", lw = 7, la = 0.7)
+plot!(nc_Y; vars = (:param, :x), label = "dY/dt = 0", lw = 7, la = 0.7)
+for ss in sss
+    color, markershape = if steady_state_stability(ss, bs_switch, ps)
+        :blue, :circle
+    else
+        :red, :star4
+    end
+    scatter!((ss[2], ss[1]); color, markershape, label = "", markersize = 10)
+end
+scatter!([], []; color = :blue, markershape = :circle, label = "Stable stead state")
+scatter!([], []; color = :red, markershape = :star4, label = "Unstable stead state")
+plot!(xlimit = span, ylimit = span, xguide = "X", yguide = "Y", legendfontsize = 10, size = (500,500))
+```
+Here we can see how the steady states occur at the nullclines intersections.
+
+## [Plotting system directions in phase space](@id nullcline_plotting_directions)
+
+```@example nullcline_plotting
+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
+function plot_xy_arrow!(Xval, Yval)
+    dX, dY = get_XY(Xval, Yval)
+    dX = dX > 0 ? 0.05 : -0.05
+    dY = dY > 0 ? 0.05 : -0.05
+    plot!([Xval, Xval + dX], [Yval, Yval]; color = :black, lw = 2, arrow = :arrow, label = "")
+    plot!([Xval, Xval], [Yval, Yval + dY]; color = :black, lw = 2, arrow = :arrow, label = "")
+end
+arrow_positions = [(0.25, 0.25), (0.75, 0.75), (0.35, 0.8), (0.8, 0.35), (0.02, 1.1), (1.1, 0.02)]
+foreach(pos -> plot_xy_arrow!(pos...), arrow_positions)
+plot!()
+```