Plotting Nullclines and Steady States in Phase Space

In this tutorial we will show how to extract a system's steady states and nullclines, and how to plot these in phase space. Generally, while nullclines are not directly needed for most types analysis, plotting these can give some understanding of a system's steady state and stability properties.

For an ordinary differential equation

\[\begin{aligned} \frac{dx_1}{dt} &= f_1(x_1, x_2, ..., x_n) \\ \frac{dx_2}{dt} &= f_2(x_1, x_2, ..., x_n) \\ &\vdots\\ \frac{dx_n}{dt} &= f_n(x_1, x_2, ..., x_n) \\ \end{aligned}\]

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

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.

using Catalyst
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

\[ \begin{align*} \varnothing &\xrightarrow{\frac{K^{n} v}{K^{n} + Y^{n}}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{K^{n} v}{K^{n} + X^{n}}} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{1.0} \varnothing \\ \mathrm{Y} &\xrightarrow{1.0} \varnothing \end{align*} \]

Next, we compute the steady states using homotopy continuation.

import HomotopyContinuation
ps = [:v => 1.0, :K => 0.6, :n => 4.0]
sss = hc_steady_states(bs_switch, ps; show_progress = false)

3-element Vector{Vector{Float64}}:
 [0.9986388662510188, 0.11528540639768478]
 [0.563112248656171, 0.563112248656171]
 [0.11528540639768477, 0.9986388662510189]

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.

nlprob = NonlinearProblem(bs_switch, [:X => 0.0, :Y => 0.0], ps)
function get_XY(Xval, Yval)
    prob = Catalyst.remake(nlprob; u0 = [:X => Xval, :Y => Yval])
    return nlprob.f(prob.u0, prob.p)
end

Next, we plot our nullclines by using Plot.jl's contour function. 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 (and use this to determine how to plot the steady state markers).

using Plots
# 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
    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")

# Finishing touches.
plot!(xguide = "X", yguide = "Y", xlimit = (0.0, 1.25), ylimit = (0.0, 1.25), legendfontsize = 10, size = (600,600), legend = :topright)

Here we can see how the steady states occur at the nullclines intersections.

Plotting system directions in phase space

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.

# 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)
    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

# Plots the ODE's direction in phase space at selected positions.
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!()

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.