Logo

Analyzing critical transitions with Python

Scheffer et al. (2001) describe the following model as a "minimal model of an ecosystem showing hysteresis":

  • $\frac{dx}{dt}=a-bx+rf(x)$ (1)

with the Hill-function (see below)

  • $f(x)=\frac{x^p}{x^p + h^p}$ (2)

This equation is used in Carpenter et al. (1999) to outline the essential dynamics of eutrophication in lakes. According to Scheffer et al. (2001), a model interpretation of equation (1) and (2) could be to regard $x$ as a nutrient suspended in phytoplankton in a lake causing turbidity. $a$ would indicate the loading rate of the nutrient, $b$ its removal rate, and $r$ could be an internal recycling rate of the nutrient (see also Critical Transitions).

In the following, this model will be discussed in the framework of an iPython-notebook that provides Python code for further experimenatation. In the beginning, the focus lies on the varition of the parameter $a$ which causes alternative equilibria and has further implications for the behavior of the system.

The Hill-function

As a special logistic function the Hill-function allows to choose the position ($h$) where the development reaches half of its saturation. The steepness of its slope is determined through the exponent $p$.

In [2]:
#the Hill function
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
fig = plt.figure(figsize=(10,4))

h = 1

for p in range(2,14,2):
    X = []; FX = []
    for x in np.linspace(0,3,100):
        fx = x**p/(x**p + h**p)
        X.append(x)
        FX.append(fx)    
    plt.plot(X, FX, label = 'p = ' + str(p))
    
plt.plot([h,h],[0,1], 'k--', label = 'h (='  + str(h) + ')')
plt.title('The Hill-function')
plt.xlabel('$x$', fontsize=18)
plt.ylabel('$f(x)$', fontsize=18)
plt.grid()
plt.text(0.1, .9, 'fig. 1', fontsize=14)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
Out[2]:
<matplotlib.legend.Legend at 0x74e6470>

The behavior of $x$ in time

In this notebook the differential equation (1) will be treated as a difference equation. Alternatively scipy ode solver or other tools could be used for numerically calculating differential equations. Transforming the differential equation (1) (with respect to (2)) into a difference equation gives

  • $x_{t+1} = x_t + \Delta x_t$ (3)

with

  • $\Delta x_t = a-bx_t+\frac{x_t^p}{x_t^p + h^p}$ (4)

Using equation (3), $x(t)$ can be calculated with a for-loop in python, simply iterating through the whole range of $t$. The parameters are set to $ $ $h$ = 1.0, $b$ = 1.0, $p$ = 12 and $r$ = 0.6, $a$ varies in the range from 0.5 to 0.9 with an increment of $10^{-4}$. The limitation for time $t$ is not that important in difference equations, but the number of calculations and their increment is. Here we use 2000 steps from 0.5 to 1. Additionally $\Delta x_t$ is set to $\frac{\Delta x_t}{100}$ to smooth the curve.

In [3]:
#the behavior of x(t)

#load necessary python-packages:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

#define figure properties for upcoming plots:
fig = plt.figure(figsize=(15,7))
mpl.rc('font', size = 14)

#define parameters and predefine array variables:
h = 1.; b = 1.; p = 12; r = 0.6; A=[]; XA=[]

#vary parameter a in a range of (0.4,0.9):
a = 0.5
while a <= 0.9:
    x = 1; X = []; T = []
    #calculate x(t) for parameter t in a range of (0.5,1):
    for t in np.linspace(0.5,1,2000):
        #the difference equation:
        delta_x = a - b * x + r *(x**p / (x**p + h**p))
        x += delta_x/100
        X.append(x); T.append(t)
    #plot x(t) for specific values of parameter a:
    if len(str(a)) == 3:
        plt.subplot(2, 1, 1)   
        plt.plot(T,X,label = 'a = ' + str(a))
    A.append(a)
    XA.append(x)
    a += 0.0001

#configure plot properties for subplot(2, 1, 1):
plt.tight_layout(h_pad=5)
plt.title('x(t) for specific values of a')
plt.axis((0.5,1,0.3,1.7))
plt.xlabel('$t$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.grid()
plt.text(0.51, 1.5, 'fig. 2', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#plot x(a) at t=1 and configure plot properties:
plt.subplot(2, 1, 2)
plt.title('equilibria of x(a) at t = 1 for x0 = 1')
plt.plot(A,XA,color = 'black', label = 't = 1', linewidth=2, linestyle = '-')
plt.xlabel('$a$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.grid()
plt.text(0.51, 1.4, 'fig. 3', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
Out[3]:
<matplotlib.legend.Legend at 0x8e64110>

The upper one of the two plots above (fig. 2) shows the temporal development of $x(t)$ for specific values of $a$. Obviously, depending on parameter $a$, different equilibria (stable states) are achieved at $t$ = 1.
The lower plot (fig. 3) illustrates the development of $x$ with different values of $a$ at $t$ = 1. At $a$ = 0.7 apparently a sudden shift happens. This shift marks the tipping point for a critical transition. Note that this applies to an initial value of $x$ = 1. Different initial values can result in different equilibria, as will be discussed below.

Equilibria (= extrema) of $x$ in time

The right side of equation (1) and respectively (2)

  • $\frac{dx}{dt} = a-bx+\frac{x^p}{x^p + h^p}$ (5)

can be assigned to a function of $x$:

  • $g(x) = a-bx+\frac{x^p}{x^p + h^p}$ (6)

which can be used to analyze the behavior of $\frac{dx}{dt}$ for different values of $x$.

$\frac{dx}{dt} = 0$ indicates zero development and thus an equilibrium or extremum of $x(t)$. The points where $g(x) = 0$ represent such extrema. Since $x$ in equation (6) cannot be expressed explicitly (at least not in a convenient way), the points with $g(x) = 0$ are calculated numerically below, with a precision of $\pm 10^{-4}$.

The derivative of $g(x)$ with respect to $x$ allows to determine whether an extremum represents a stable or an unstable state of the system:

  • $\frac{dg}{dx} < 0$ ... stable state
  • $\frac{dg}{dx} > 0$ ... unstable state
  • $\frac{dg}{dx} = 0$ ... higher derivatives must be taken into account
In [4]:
#equilibria of x(t)

fig = plt.figure(figsize=(15,7))
mpl.rc('font', size = 14)

EX1 = []; EX2 = []; EX3 = []

a = 0.5
while a <= 0.9:
    GX = []; X = []
    #calculate g(x) for parameter t in a range of (0,2):
    for x in np.linspace(0,2,1000):
        gx = a-b*x + (r*x**p)/(x**p + h**p)
        GX.append(gx); X.append(x)
        if -10**-4 < gx < 10**-4:
            #calculate g'(x) (= dg(x)/dx) to "check" stability:
            Dgx = -b+(r*p*x**(p-1)*h**p)/((x**p + h**p)**2)
            #write x-coordinate of extremum and depending value of parameter a in an array:
            if Dgx > 0: EX2.append([x,a])      #array for unstable extrema
            if Dgx < 0:                        
                if x < 1: EX1.append([x,a])    #array for stable extrema with x<1
                if x > 1: EX3.append([x,a])    #array for stable extrema with x>1
    #plot f(x) for specific values of a:
    if len(str(a)) == 3:        
        plt.subplot(2, 1, 1)
        plt.plot(X,GX,label = 'a = ' + str(a))
    a += 0.0001

#configure plot properties for subplot(2, 1, 1):
plt.subplot(2, 1, 1);
plt.tight_layout(h_pad=5)
plt.title('g(x) for specific values of a')
plt.xlabel('$x$', fontsize=18)
plt.ylabel('$g_{(x)}$', fontsize=18)
#plt.plot(0, 'k')
plt.grid()
plt.text(0.05, -0.9, 'fig. 4', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#sort equilibria containing arrays in respect to the increasing order of parameter a:
EX1 = zip(*sorted(EX1,key=lambda x: x[0])); EX1x = EX1[1]; EX1y = EX1[0]
EX2 = zip(*sorted(EX2,key=lambda x: x[0])); EX2x = EX2[1]; EX2y = EX2[0]
EX3 = zip(*sorted(EX3,key=lambda x: x[0])); EX3x = EX3[1]; EX3y = EX3[0]

#plot curve representing equilibria of x depending on the variation of parameter a:
plt.subplot(2, 1, 2);
#lower part of the curve (stable extrema):
plt.plot(EX1x,EX1y,color = 'k', linewidth=2, label = 'stable equilibria');
#central part of the curve (unstable extrema):
plt.plot(EX2x,EX2y,color = 'k', linewidth=2, label = 'unstable equilibria', linestyle = 'dashdot');
#upper part of the curve (stable extrema):
plt.plot(EX3x,EX3y,color = 'k', linewidth=2);

#configure plot properties:
plt.title('equilibria of x(t) for specific values of a')
plt.grid();
plt.xlabel('$a$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.text(0.51, 1.4, 'fig. 5', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
268
Out[4]:
<matplotlib.legend.Legend at 0x92c2530>

Figure 4 shows an array of curves indicating $g(x)$ for specific values of $a$. The $x$-values of the points, where these curves cross the zero-line, indicate extrema of $x(t)$.

Figure 5 shows the $x$-coordinates of the $g(x)=0$-points in respect to $a$, calculated for $a$ in a range from 0.5 to 0.9 with an increment of $10^{-4}$.
The calculation of $g'(x)$ $(=\frac{dg}{dx})$ provides information about the stability of an extremum, indicated by dashed and normal linestyles.

Hysteresis and critical transition

In [2]:
#hysteresis plotting function

import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline

def hysteresis(EX1x,EX1y,EX2x,EX2y,EX3x,EX3y,x_label):

    fig = plt.figure(figsize=(15,7))
    mpl.rc('font', size = 14)
    ax = fig.add_subplot(1,1,1)
    
    #plot curve containing the extrema, splitted in lower (solid line), central (dashdotted line)
    #and upper part (solid line) in black color:
    ax.plot(EX1x,EX1y,color = 'k', linewidth=3, label = 'stable equilibria')
    ax.plot(EX2x,EX2y,color = 'k', linewidth=3, linestyle = 'dashdot', label = 'unstable equilibria')
    ax.plot(EX3x,EX3y,color = 'k', linewidth=3)
    
    #plot development of x_extremum starting at maximum value of parameter a
    #find x-coordinate in EX3x at which critical transition occurs:
    for index3, value3 in enumerate(EX3x):
        if -0.005 < (EX1x[-1] - value3) < 0.005:
            break
    line, = ax.plot(EX1x,EX1y,color = 'c', linewidth=8, linestyle='-',
            label='development of $x_{equilibrium}$ for increasing \nvalues of '
            + x_label + ', starting at ' + x_label + '$_{min}$')
    #line.set_dashes((dash_length,dash_distance)) allows to vary the dash-style
    line.set_dashes((10,20))
    line, = ax.plot([EX1x[-1],EX1x[-1]],[EX1y[-1],EX3y[index3]],color = 'c', linewidth=8, linestyle='-')
    line.set_dashes((10,20))
    line, = ax.plot(EX3x[index3:],EX3y[index3:],color = 'c', linewidth=8, linestyle='-')
    line.set_dashes((10,20))
    
    #plot development of x_extremum starting at minimum value of parameter a
    #find x-coordinate in EX1x at which critical transition occurs:
    for index1, value1 in enumerate(EX1x):
        if -0.005 < (EX3x[0] - value1) < 0.005:
            break
    line, = ax.plot(EX3x,EX3y,color = 'r', linewidth=8, linestyle='-',
            label='development of $x_{equilibrium}$ for decreasing \nvalues of '
            + x_label + ', starting at ' + x_label + '$_{max}$')
    line.set_dashes((5,20))
    line, = ax.plot([EX3x[0],EX3x[0]],[EX3y[0],EX1y[index1]],color = 'r', linewidth=8, linestyle='-')
    line.set_dashes((5,20))
    line, = ax.plot(EX1x[:index1],EX1y[:index1],color = 'r', linewidth=8, linestyle='-')
    line.set_dashes((5,20))

    return ax, index3, index1
In [5]:
#plot hysteresis

ax, index3, index1 = hysteresis(EX1x,EX1y,EX2x,EX2y,EX3x,EX3y,'$a$')

#add arrows to the graph:
#cyan arrows:
ax.annotate('', xy=(0.85, 1.35), xytext=(0.8, 1.3),arrowprops=dict(facecolor='c',lw=0),)
ax.annotate('', xy=(0.75, 0.7), xytext=(0.7, 0.63),arrowprops=dict(facecolor='c',lw=0),)
ax.annotate('', xy=(0.58, 0.5), xytext=(0.53, 0.45),arrowprops=dict(facecolor='c',lw=0),)
#red arrows:
ax.annotate('', xy=(0.8, 1.46), xytext=(0.85, 1.52),arrowprops=dict(facecolor='r',lw=0),)
ax.annotate('', xy=(0.67, 1.3), xytext=(0.72, 1.37),arrowprops=dict(facecolor='r',lw=0),)
ax.annotate('', xy=(0.53, 0.61), xytext=(0.58, 0.66),arrowprops=dict(facecolor='r',lw=0),)
#black arrows: 
ax.annotate('', xy=(0.63, 0.65), xytext=(0.63, 1.6), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.annotate('', xy=(0.66, 1.22), xytext=(0.66, 1.6), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.annotate('', xy=(0.63, 0.62), xytext=(0.63, 0.45), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.annotate('', xy=(0.66, 1.18), xytext=(0.66, 1.07), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.annotate('', xy=(0.66, 0.68), xytext=(0.66, 1.05), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.annotate('', xy=(0.66, 0.64), xytext=(0.66, 0.45), arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)

#configure plot properties:
ax.set_title('Hysteresis depicting the development of equilibria (x-cor) depending on variation of a')
ax.grid();
ax.set_xlabel('$a$',fontsize=18);
ax.set_ylabel('$x$',fontsize=18);
ax.set_xlim([EX1x[0],EX3x[-1]])
ax.set_ylim([EX1y[0]-0.1,EX3y[-1]+0.2])
ax.text(0.51, 1.6, 'fig. 6', fontsize=18)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
Out[5]:
<matplotlib.legend.Legend at 0xbc0bac8>
In [5]:
#tipping points

print 'tipping point for increasing parameter a:'
print 'a =', round(EX1x[-1],2), ', x_low =', round(EX1y[-1],2), ', x_high =', round(EX3y[index3],2)
print 'tipping point for decreasing parameter a:'
print 'a =', round(EX3x[0],2), ', x_low =', round(EX1y[index1],2), ', x_high =', round(EX3y[0],2)
tipping point for increasing parameter a:
a = 0.78 , x_low = 0.86 , x_high = 1.36
tipping point for decreasing parameter a:
a = 0.64 , x_low = 0.64 , x_high = 1.13

Figure 6 highlights the development of equilibria in the system. The difference between blue development with increasing $a$-values and red development with decreasing $a$-values indicates the hysteresis in the system:

The dotted blue curve illustrates the development of the equilibrium starting from $a$-values $\lesssim$ 0.64 and increasing $a$ steadily. For $a \lesssim$ 0.78 $x$ behaves predictable. At $a \sim$ 0.78 the tipping point is reached and the development shifts from $x$ = 0.86 to $x$ = 1.36. This shift represents a critical transition, for example from a low to a high nutrient load ($=x$) in a lake. Further increasing $a$ lets $x$ behave predictable.
Decreasing $a$, now that it is $\gtrsim$ 0.78, does not result in the same transition. The development of $x$ for decreasing values of $a$ is illustrated with the dash-dotted red curve. Down to a value of $a$ = 0.64, $x$ decreases steadily. At $a$ = 0.64 the second tipping point is achieved and a shift from $x$ = 1.13 to $x$ = 0.64 occurs. For $a \lesssim$ 0.64 $x$ behaves predictable.

This effect is called hysteresis: The system behaves in respect to its history. The vertical thin black arrows In figure 6 indicate the attraction of equilibria in neighboring regions around the tipping point at $a \sim$ 0.64. Examples for this difference in the attraction behavior are given in the next paragraph.

Different attraction behavior and the meaning of stability

In [6]:
#different attraction behavior

fig = plt.figure(figsize=(15,10))
mpl.rc('font', size = 14)

#plot the curves representing equilibria of x depending on the variation of parameter a:
ax = fig.add_subplot(2,1,1)
ax.plot(EX1x,EX1y,color = 'k', linewidth=2, label = 'stable equilibria')
ax.plot(EX2x,EX2y,color = 'k', linewidth=2, label = 'unstable equilibria', linestyle = 'dashdot')
ax.plot(EX3x,EX3y,color = 'k', linewidth=2)

ax.plot(0.7,1.05,color = 'g',marker='o')
ax.annotate('', xy=(0.7, 1.25), xytext=(0.7, 1.07),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.plot(0.7,1.27,color = 'g',marker='x',ms=8,mew=3)

ax.plot(0.7,1,color = 'c',marker='o')
ax.plot(0.7,1,color = 'c',marker='x',ms=8,mew=3)

ax.plot(0.7,0.95,color = 'm',marker='o')
ax.annotate('', xy=(0.7, 0.73), xytext=(0.7, 0.92),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.plot(0.7,0.72,color = 'm',marker='x',ms=8,mew=3)

ax.plot(0.775,0.95,color = 'r',marker='o')
ax.annotate('', xy=(0.775, 1.34), xytext=(0.775, 0.97),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.plot(0.775,1.36,color = 'r',marker='x',ms=8,mew=3)

ax.plot(0.775,0.5,color = 'b',marker='o')
ax.annotate('', xy=(0.775, 0.81), xytext=(0.775, 0.53),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.plot(0.775,0.85,color = 'b',marker='x',ms=8,mew=3)

ax.plot(0.785,0.5,color = 'y',marker='o')
ax.annotate('', xy=(0.785, 1.33), xytext=(0.785, 0.53),arrowprops=dict(arrowstyle="->",connectionstyle="arc3",lw=1),)
ax.plot(0.785,1.37,color = 'y',marker='x',ms=8,mew=3)

#configure plot properties:
ax.set_title('Achieved equilibria of x for specific values of x0 and a')
ax.grid();
ax.set_xlabel('$a$', fontsize=18)
ax.set_ylabel('$x$', fontsize=18)
ax.text(0.51, 1.4, 'fig. 7', fontsize=18)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#define parameters and predefine array variables:
h = 1.; b = 1.; p = 12; r = 0.6; A=[]; XA=[]
start_xy_list = [[0.7,1.05,'g'],[0.7,1.0,'c'],[0.7,0.95,'m'],[0.775,0.95,'r'],[0.775,0.5,'b'],[0.785,0.5,'y']]

#vary parameter a in a range of (0.4,0.9):
for start_xy in start_xy_list:
    a = start_xy[0]
    x = start_xy[1]
    color = start_xy[2]
    X = []; T = []
    for t in np.linspace(0.5,0.75,2500):
        #the difference equation:
        delta_x = a - b * x + r *(x**p / (x**p + h**p))
        x += delta_x/100
        X.append(x); T.append(t)
    plt.subplot(2, 1, 2)   
    plt.plot(T,X,label = '$P$ $(a,x_0) = $'+ str(start_xy[:2]),color = color, linewidth=2)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

#configure plot properties:
plt.tight_layout(h_pad=3)
plt.title('Development of x(t) for specific values of x0 and a (time domain)')
plt.grid();
plt.xlim(0.5,0.75)
plt.ylim(0.5,1.5)
plt.xlabel('$t$', fontsize=18)
plt.ylabel('$x$', fontsize=18)
plt.text(0.51, 1.3, 'fig. 8', fontsize=18)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
Out[6]:
<matplotlib.legend.Legend at 0xbf942b0>

Figure 7 illustrates the attraction conditions of $x$ for selected initial values ($x_0$) and different values of $a$. The initial values of $x$ ($a$-values stay constant) are represented by circles. The crosses indicate the final values to which these dynamics are attracted (equilibria, stable states). Black arrows indicate the development over time. The list of initial conditions in the legend to figure 8 as well as the final values of $x$ apply to both plots.

Figure 8 illustrates the development of $x$ in time. The value of x for the cyan curve remains constant due to the fact that its intial value is set on an equilibrium. As can be seen in fig. 7, the point $P($0.7,1$)$ ($a$ = 0.7, $x$ = 1) represented by the cyan circle lais on the unstable equilibria of the system. It is a mathematical particularity. In nature, even the smallest deviation would suffice to detract the dynamics to one of the stable equilibria.
This can be seen at the magenta circle in fig. 7, which is positioned only slightly below the cyan circle. At the point $P($0.7,0.95$)$ however, $x$ gets attracted by the lower stable equlibrium.

A partly counterintuitive and therefore very interesting behavior shows the yellow curve in figure 8. In fig. 7, the development of the the yellow circle at $P($0.785,0.5$)$ is straightforward. In figure 8 it is not so clear. At $t \sim$ 0.54 the yellow curve flattens first, but then shows more or less constant growth of $x$ up to a value of $t$ = 0.6. At $t \sim$ 0.6 a tipping point is achieved and $x$ starts to increase again, before seemingly logistically approaching a stable state at $t$ > 0.7. This intriguing behavior of the development of $x$ for the initial condition $P($0.785,0.5$)$ can be understood when regarding the attraction conditions as a hypersurface, which is shown in the next plot below.

Illustration of attraction conditions as hypersurface

In [11]:
#Illustration of attraction conditions as hypersurface

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.tri as mtri

#g(x)
X = np.linspace(0.5,1.5,25)
A = np.linspace(0.5,0.9,25)
X, A = np.meshgrid(X, A)
GX = np.array(A-b*X + (r*X**p)/(X**p + h**p))

#curve representing equilibria
X1 = np.linspace(0.5,EX1y[-1],100)
X2 = np.linspace(EX1y[-1],EX3y[0],100)
X3 = np.linspace(EX3y[0],1.5,100)
A1 = b*X1 - (r*X1**p)/(X1**p + h**p)
A2 = b*X2 - (r*X2**p)/(X2**p + h**p)
A3 = b*X3 - (r*X3**p)/(X3**p + h**p)

#yellow curve
Xyc = np.linspace(0.5,1.37,25)
Ayc = np.ones(len(Xyc))*0.785
yellow_curve = np.array(Ayc-b*Xyc + (r*Xyc**p)/(Xyc**p + h**p))


#plot
fig = plt.figure(figsize = (15, 10))
ax = fig.add_subplot(111, projection='3d')

ax.plot_trisurf(X.ravel(),A.ravel(),GX.ravel(), lw=0.2, cmap=cm.seismic)
m = cm.ScalarMappable(cmap=cm.seismic)
m.set_array(GX)
plt.colorbar(m, shrink=0.5, aspect=8)

ax.plot(X1,A1,0*X1,color='k',lw=1.5)
ax.plot(X2,A2,0*X2,color='k',lw=1.5,ls='--')
ax.plot(X3,A3,0*X3,color='k',lw=1.5)

ax.plot(Xyc,Ayc,yellow_curve,color='yellow',lw=2.5,ls='-')

ax.set_xticks([0.5,0.75,1,1.25,1.5])
ax.set_yticks([0.5,0.6,0.7,0.8,0.9])
ax.set_zticks([-0.5, -0.25, 0, 0.25, 0.5])
ax.invert_yaxis()

ax.set_xlabel('x', fontsize = 15)
ax.set_ylabel('a', fontsize = 15)
ax.set_zlabel('g(a,x)', fontsize = 15)
ax.text(1, 0, 0, 'fig. 9', fontsize=18)

ax.view_init(30, 120)

The hypersurface in fig. 9 illustrates $g(a,x)$ for all combinations of $a$ in range (0.5,0.9) and $x$ in range (0.5,1.5). The black curve with the dashed part in the middle represents the equilibria of the system. The red part of the hypersurface indicates the region with $g > 0$ (the darker the higher the value of $g$). The blue region represents $g < 0$ (the darker the lower the value of $g$).
As explained above, $g \simeq \frac{dx}{dt}$, therefore $g$ is the gradient of $x$. As a result, the red region indicates that if $x$ is located there it will increase. If located in the blue region, $x$ decreases, both until a stable equilibrium is acchieved.

The yellow curve in fig. 9 has the same initial condition at $P($0.785,0.5$)$ as in fig. 7 and 8. Since it starts out with a high value of $g$ (fig. 9), it increases in time (fig. 8). However, when $g$ decreases almost to zero at $x \sim 0.85$, the increase is slow, as shown in fig. 8, - and so on.
In this way, figure 9 can explain how equilibria are achieved for different initial conditions.

Hysteresis of equilibria in respect to parameters $b$ and $r$

Hysteresis can also be observed in respect to other parameters in this model, i.e. in respect to $b$ for the removal rate of the nutrient and in respect to $r$ for its recycling rate. This is shown in the following four notebook parts.

In [3]:
#hysteresis of equilibria - parameter b

import numpy as np

h = 1.; b = 0.5; p = 12; r = 0.7; a = 0.7
EX1 = []; EX2 = []; EX3 = []

while b <= 1.5:
    X = []; T = []
    for x in np.linspace(0.5,2,3000):
        fx = a-b*x + (r*x**p)/(x**p + h**p)
        X.append(fx); T.append(x)
        if -0.0001 < fx < 0.0001:
            fx_strich = -b+(r*p*x**(p-1)*h**p)/((x**p + h**p)**2)
            if fx_strich > 0: EX2.append([x,b])
            if fx_strich < 0:
                if x < 1: EX1.append([x,b])          
                if x > 1: EX3.append([x,b])             
    b += 0.0005

EX1 = zip(*sorted(EX1,key=lambda x: x[0])); EX1x = EX1[1]; EX1y = EX1[0]
EX2 = zip(*sorted(EX2,key=lambda x: x[0])); EX2x = EX2[1]; EX2y = EX2[0]
EX3 = zip(*sorted(EX3,key=lambda x: x[0])); EX3x = EX3[1]; EX3y = EX3[0]

ax, index3, index1 = hysteresis(EX1x,EX1y,EX2x,EX2y,EX3x,EX3y,'$b$')
ax.annotate('', xy=(1.0, 1.52), xytext=(0.95, 1.6),arrowprops=dict(facecolor='r',lw=0),)
ax.annotate('', xy=(1.0, 0.6), xytext=(1.05, 0.56),arrowprops=dict(facecolor='c',lw=0),)
ax.set_title('Hysteresis depicting the development of equilibria (x-cor) depending on variation of b')
ax.grid();
ax.set_xlabel('$b$', fontsize=18);
ax.set_ylabel('$x$', fontsize=18);
ax.set_xlim([EX3x[-1],EX1x[0]])
ax.set_ylim([EX1y[0]-0.1,EX3y[-1]+0.2])
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
Out[3]:
<matplotlib.legend.Legend at 0x7b51510>