Chapter 6

Nonlinear neurodynamics and bifurcations

Code by : Abolfazl Ziaeemehr - https://github.com/Ziaeemehr

# uncomment and run this line to install the package on colab
# !pip install "git+https://github.com/Ziaeemehr/spikes.git" -q

import sympy
import sympy as sp
import numpy as np
import matplotlib as mpl
from scipy.linalg import eig
import matplotlib.pyplot as plt
from IPython.display import display, Math
from spikes.solver import solve_system_of_equations
import warnings
warnings.filterwarnings("ignore")

sympy.init_printing()

LABELSIZE = 13
plt.rc('axes', labelsize=LABELSIZE)
plt.rc('axes', titlesize=LABELSIZE)
plt.rc('figure', titlesize=LABELSIZE)
plt.rc('legend', fontsize=LABELSIZE)
plt.rc('xtick', labelsize=LABELSIZE)
plt.rc('ytick', labelsize=LABELSIZE)
mpl.rcParams['image.cmap'] = 'jet'

To plot the nullclines of a two-dimensional system of differential equations, you can use Python along with libraries like NumPy, Matplotlib, and SymPy. Here’s a brief overview of how to do it:

Step 1: Define the System of Differential Equations

A two-dimensional system of differential equations can be written as:

\[\begin{split}\begin{cases} \frac{dx}{dt} = f(x, y) \\ \frac{dy}{dt} = g(x, y) \end{cases}\end{split}\]

Step 2: Find the Nullclines

• Nullcline for \(\frac{dx}{dt} = 0\) : Set \(f(x, y) = 0\) .

• Nullcline for \(\frac{dy}{dt} = 0\) : Set \(g(x, y) = 0\) .

Step 3: Use SymPy to Solve for Nullclines

You can use SymPy to symbolically solve for \(y\) in terms of \(x\) (or vice versa) where these equations are zero.

Step 4: Plot the Nullclines Using Matplotlib

After finding the expressions for the nullclines, you can use Matplotlib to plot them over a grid of values for \(x\) and \(y\).

Example

Consider the following system of differential equations:

\[\begin{split}\begin{cases} \frac{dx}{dt} = x^2 - y \\ \frac{dy}{dt} = y - x \end{cases}\end{split}\]

Explanation

1. Defining the system: The functions \(f(x, y)\) and \(g(x, y)\) represent the right-hand side of the system of differential equations.

2. Finding nullclines: We solve \(f(x, y) = 0\) and \(g(x, y) = 0\) for \(y\).

3. Plotting: The nullclines are plotted using Matplotlib, with one set for \(\frac{dx}{dt} = 0\) and another for \(\frac{dy}{dt} = 0\).

Tools Used

• SymPy: For symbolic computation, solving the equations for nullclines.

• NumPy: For numerical operations and creating arrays of values.

• Matplotlib: For plotting the nullclines on a 2D plane.

This approach helps visualize where the system’s derivatives are zero, giving insight into the system’s behavior, such as identifying equilibrium points and possible phase portrait structures.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
from ipywidgets import interact
from scipy.integrate import odeint
from spikes.plot import plot_nullclines

# Example usage:
x, y = sp.symbols('x y')
f = x**2 - y
g = y - x

fig, ax = plt.subplots(figsize=(4,4))
plot_nullclines(f, g, symbols=['x', 'y'], ax=ax);

# adding direction fields
def fp(x, y):
    return x**2 - y

def gp(x, y):
    return y - x


def plot_vector_field(
    f, g, ax,
    xrange=(-2, 2),
    yrange=(-2, 2),
    num_points=20,
    rotatexy=False):
    x = np.linspace(xrange[0], xrange[1], num_points)
    y = np.linspace(yrange[0], yrange[1], num_points)
    X, Y = np.meshgrid(x, y)

    u = f(X, Y)
    v = g(X, Y)

    norm = np.sqrt(u**2 + v**2)
    u_norm = np.real(u / norm)
    v_norm = np.real(v / norm)
    if rotatexy:
        ax.quiver(Y, X, v_norm, u_norm);
        # ax.streamplot(Y, X, v_norm, u_norm, color='k', linewidth=0.5, cmap=plt.cm.autumn)
    else:
        ax.quiver(X, Y, u_norm, v_norm)
        # ax.streamplot(X, Y, u_norm, v_norm, color='k', linewidth=0.5, cmap=plt.cm.autumn)

x, y = sp.symbols('x y')
f = x**2 - y
g = y - x

fig, ax = plt.subplots(figsize=(4,4))
plot_nullclines(f, g, symbols=['x', 'y'], ax=ax);
plot_vector_field(fp, gp, ax);

adding an interactive trajectory solution

If you are getting messy output just comment the line iteract()

def ode(x, t):
    return np.array([fp(x[0], x[1]), gp(x[0], x[1])])

def plot_trajectory(x0=-1., y0=0, tspan=2):

    fig, ax = plt.subplots(figsize=(5, 5));
    plot_nullclines(f, g, symbols=['x', 'y'], ax=ax);
    plot_vector_field(fp, gp, ax)
    t_range = np.arange(0, tspan, 0.1)

    sol = odeint(ode, y0=np.array([x0,y0]), t=t_range)
    ax.plot(sol[:,0], sol[:, 1] , c='g', alpha=.5);
    ax.plot([x0], [y0], 'go', ms=10, alpha=0.5);
    ax.set_title(f'Phase Portrait (x0={x0:.2f}, y0={y0:.2f})')
    # plt.show()

# comment the following line if you don't want to use interactive mode:
# interact(plot_trajectory, x0=(-1.99,1.99,.1), y0=(-1.99,1.99,.1), tspan=(1,5,0.5))

# uncomment the following for using without interactive mode
plot_trajectory()

Equation 6.3

\[\begin{split}\begin{align*} \frac{dA}{dt} &= \frac{1}{\tau_A}\big( -A + 2B \big) \\ \frac{dB}{dt} &= \frac{1}{\tau_B}\big( -B + \frac{L}{1+A} \big) \\ \end{align*}\end{split}\]

def fp(A, B):
    return -A + 2 * B

def gp(A, B):
    return -B + L / (A+1)

def ode(x, t):
    return np.array([fp(x[0], x[1]), gp(x[0], x[1])])

A, B = sp.symbols('A B')
L = 10
f = -A + 2 * B
g = -B + L / (A+1)

def plot_trajectory(A0=4, B0=6, tspan=3):

    fig, ax = plt.subplots(figsize=(5, 5));
    plot_nullclines(g, f, symbols=['B', 'A'], x_range=(0,10), y_range=(0,10), ax=ax);
    plot_vector_field(fp, gp, ax, xrange=(0,10), yrange=(0,10), rotatexy=True)
    t_range = np.arange(0, tspan, 0.1)

    sol = odeint(ode, y0=np.array([A0,B0]), t=t_range)
    ax.plot(sol[:,1], sol[:, 0] , c='k', alpha=.5);
    ax.plot([B0], [A0], 'ko', ms=10, alpha=0.5);
    ax.set_title(f'Phase Portrait (A0={A0:.2f}, B0={B0:.2f})')


# comment the following line if you don't want to use interactive mode:
# interact(plot_trajectory, A0=(0.1,9.99,.1), B0=(0.1,9.99,.1), tspan=(1,10,0.5))

# uncomment the following for using without interactive mode
plot_trajectory()

Related tool

pplane

There is also a java/matlab tool available for visualization of the nullclines and vector fields here.

Stabiliby of steady states

The process includes the calculation of the Jacobian matrix, finding the eigenvalues of the Jacobian matrix at the equilibrium point, and classifying the equilibrium point based on the nature of the eigenvalues.

We’ll use the sympy library for symbolic differentiation (to compute the Jacobian) and numpy for numerical linear algebra (to find eigenvalues).

Here’s a step-by-step Python code that implements the process:

Explanation:

1. Step 1 (Equilibrium Points):

   • We define a system of equations \(\frac{dX}{dt} = F(X)\).

   • The equilibrium points are found by solving \(F(X) = 0\).

2. Step 2 (Jacobian Matrix):

   • The Jacobian matrix is calculated symbolically using sympy.jacobian.

3. Step 3 (Eigenvalue Analysis and Classification):

   • The Jacobian matrix at each equilibrium point is substituted into and evaluated numerically.

   • The eigenvalues are calculated using numpy.linalg.eigvals.

   • The type of equilibrium point is classified based on the eigenvalues:

     • Node (Stable/Unstable): All real and either positive (unstable) or negative (stable).

     • Saddle Point: Mixed signs among real eigenvalues.

     • Spiral (Stable/Unstable): Complex eigenvalues with negative or positive real parts.

     • Center: Purely imaginary eigenvalues (neutral stability).

Customizing the Code:

1. System of Equations: Modify F for your system.

2. Variables: If your system has more variables, extend the variables x1, x2, etc.

3. Jacobian Calculation: The Jacobian calculation automatically scales with the number of variables in the system.

import sympy as sp
import numpy as np
from IPython.display import display, Math
sp.init_printing()

def find_jacobian(F, X):
    """
    Compute the Jacobian matrix for a system of equations F(X).

    Parameters
    ----------
    F: sympy.Matrix
        The system of equations represented as a sympy Matrix.
    X: sympy.Symbol or iterable of sympy.Symbol
        The independent variables of the system of equations.

    Returns
    -------
    J: sympy.Matrix
        The Jacobian matrix of the system of equations.

    """
    J = F.jacobian(X)
    return J

def equilibrium_points(F, X):
    """
    Solve for the equilibrium points by setting F(X) = 0.

    Parameters
    ----------
    F: sympy.Matrix
        The system of equations represented as a sympy Matrix.
    X: sympy.Symbol or iterable of sympy
        The independent variables of the system of equations.

    Returns
    -------
    eq_points: sympy.Eq
        List of equilibrium points in the system.

    """
    eq_points = sp.solve(F, X)
    return eq_points

def has_complex_elements(matrix:sympy.Matrix):
    return matrix.has(sympy.I)


def classify_equilibrium(Jacobian, X_eq, X):
    """
    Classify the equilibrium point by calculating eigenvalues of the Jacobian at X_eq.

    Parameters
    ----------
    Jacobian: sympy.Matrix
        The Jacobian matrix of the system of equations.
    X_eq: sympy.Matrix
        The equilibrium point.
    X: sympy.Symbol or iterable of sympy
        The independent variables of the system of equations.

    Returns
    -------
    classification: str
        The classification of the equilibrium point.
    """
    # Substitute the equilibrium point into the Jacobian
    J_at_eq = Jacobian.subs([(x, x_eq) for x, x_eq in zip(X, X_eq)])
    # Evaluate symbolic Jacobian numerically
    J_at_eq = J_at_eq.evalf()
    # check if there is any complex value in the Jacobian
    if has_complex_elements(J_at_eq):
        J_at_eq_np = np.array(J_at_eq.tolist(), dtype=complex)
    else:
        J_at_eq_np = np.array(J_at_eq.tolist(), dtype=float)

    # Compute eigenvalues
    # J_at_eq_np = np.array(J_at_eq).astype(np.float64)  # Convert to numpy array
    eigenvalues = np.linalg.eigvals(J_at_eq_np)  # Compute eigenvalues

    # Classify based on eigenvalues
    real_parts = [ev.real for ev in eigenvalues]
    imag_parts = [ev.imag for ev in eigenvalues]

    if all(np.isreal(eigenvalues)):  # If all eigenvalues are real
        if all(ev > 0 for ev in real_parts):
            stability = "Unstable Node"
        elif all(ev < 0 for ev in real_parts):
            stability = "Stable Node"
        else:
            stability = "Saddle Point"
    else:  # Complex eigenvalues
        if all(ev < 0 for ev in real_parts):
            stability = "Stable Spiral"
        elif all(ev > 0 for ev in real_parts):
            stability = "Unstable Spiral"
        else:
            stability = "Center (Neutral Stability)"

    return {"eigenvalues": eigenvalues, "stability": stability}

# Define the variables (states)
x, y = sp.symbols('x y')
f = x**2 - y
g = y - x
F = sp.Matrix([f, g])
X = sp.Matrix([x, y])

# Step 1: Find equilibrium points
eq_points = equilibrium_points(F, X)

# Step 2: Calculate the Jacobian matrix
Jacobian = find_jacobian(F, X)
display("Jacobian:", Jacobian)

# Step 3: Classify each equilibrium point
for X_eq in eq_points:
    # display(Math(f'x_{"{eq}"} = {sp.latex(X_eq)}'))
    D = classify_equilibrium(Jacobian, X_eq, X)
    classification = D['stability']
    eigenvalues = D['eigenvalues']
    print(f"Equilibrium point {X_eq} is classified as: {classification:>20s}, eigenvalues: {eigenvalues}")

'Jacobian:'

$\displaystyle \left[\begin{matrix}2 x & -1\\-1 & 1\end{matrix}\right]$

Equilibrium point (0, 0) is classified as:         Saddle Point, eigenvalues: [-0.61803399  1.61803399]
Equilibrium point (1, 1) is classified as:        Unstable Node, eigenvalues: [2.61803399 0.38196601]

A, B = sp.symbols('A B')

L = 10
f = -A + 2 * B
g = -B + L / (A+1)
F = sp.Matrix([f, g])
X = sp.Matrix([A, B])

eq_points = equilibrium_points(F, X)
display("Equilibrium point:", eq_points)
Jacobian = find_jacobian(F, X)

display("Jacobian:", Jacobian)
for X_eq in eq_points:
    D = classify_equilibrium(Jacobian, X_eq, X)
    classification = D['stability']
    eigenvalues = D['eigenvalues']
    print(f"{'Equilibrium point: ':<30s} {X_eq}")
    print(f"{'Eigenvalues: ':<30s} {eigenvalues}")
    print(f"{'Classification: ':<30s} {classification}")
    print("-"*50)

'Equilibrium point:'

$\displaystyle \left[ \left( -5, \ - \frac{5}{2}\right), \ \left( 4, \ 2\right)\right]$

'Jacobian:'

$\displaystyle \left[\begin{matrix}-1 & 2\\- \frac{10}{\left(A + 1\right)^{2}} & -1\end{matrix}\right]$

Equilibrium point:             (-5, -5/2)
Eigenvalues:                   [-1.+1.11803399j -1.-1.11803399j]
Classification:                Stable Spiral
--------------------------------------------------
Equilibrium point:             (4, 2)
Eigenvalues:                   [-1.+0.89442719j -1.-0.89442719j]
Classification:                Stable Spiral
--------------------------------------------------

A short-term memory circuit

\[\begin{split}\begin{align*} \frac{dE_1}{dt} &= \frac{1}{\tau} \big( -E_1 + \frac{100(3E_2)^2}{120^2 + (3E_2)^2} \big) \\ \frac{dE_2}{dt} &= \frac{1}{\tau} \big( -E_2 + \frac{100(3E_1)^2}{120^2 + (3E_1)^2} \big) \end{align*}\end{split}\]

where \(\tau=20\) ms.

• find fixed points,

• plot isoclines,

• find stability characteristics of the fixed points,

tau = 20
E1, E2 = sp.symbols('E1 E2')
X = sp.Matrix([E1, E2])
F1 = 1/tau *(-E1 + (100 * (3*E2)**2)/ (120**2 + (3*E2)**2))
F2 = 1/tau * (-E2 + (100 * (3*E1)**2)/ (120**2 + (3*E1)**2))
F = sp.Matrix([F1, F2])

fig, ax = plt.subplots(1, figsize=(4,4))
plot_nullclines(F1, F2, symbols=['E1', 'E2'],
                x_range=(-10,100), y_range=(-10,100),
                num_points=100, ax=ax, loc='best')

eq_points = equilibrium_points(F, X)

Jacobian = find_jacobian(F, X)

display("Jacobian:", Jacobian)
for X_eq in eq_points:
    D = classify_equilibrium(Jacobian, X_eq, X)
    classification = D['stability']
    eigenvalues = D['eigenvalues']
    print(f"{'Equilibrium point: ':<30s} {X_eq}")
    print(f"{'Eigenvalues: ':<30s} {eigenvalues}")
    print(f"{'Classification: ':<30s} {classification}")
    print("-"*50)

'Jacobian:'

$\displaystyle \left[\begin{matrix}-0.05 & - \frac{810.0 E_{2}^{3}}{\left(9 E_{2}^{2} + 14400\right)^{2}} + \frac{90.0 E_{2}}{9 E_{2}^{2} + 14400}\\- \frac{810.0 E_{1}^{3}}{\left(9 E_{1}^{2} + 14400\right)^{2}} + \frac{90.0 E_{1}}{9 E_{1}^{2} + 14400} & -0.05\end{matrix}\right]$

Equilibrium point:             (0.0, 0.0)
Eigenvalues:                   [-0.05 -0.05]
Classification:                Stable Node
--------------------------------------------------
Equilibrium point:             (20.0000000000000, 20.0000000000000)
Eigenvalues:                   [ 0.03 -0.13]
Classification:                Saddle Point
--------------------------------------------------
Equilibrium point:             (80.0000000000000, 80.0000000000000)
Eigenvalues:                   [-0.03 -0.07]
Classification:                Stable Node
--------------------------------------------------
Equilibrium point:             (-6.89655172413793 - 13.1577820885096*I, -6.89655172413793 + 13.1577820885096*I)
Eigenvalues:                   [ 0.0577033+0.j -0.1577033+0.j]
Classification:                Saddle Point
--------------------------------------------------
Equilibrium point:             (-6.89655172413793 + 13.1577820885096*I, -6.89655172413793 - 13.1577820885096*I)
Eigenvalues:                   [ 0.0577033+0.j -0.1577033+0.j]
Classification:                Saddle Point
--------------------------------------------------

Hysteresis, bifurcation, and memory

• TODO!

Adaptation, forgetting, and catastrophe theory

\[\begin{split}\begin{align*} \frac{dE_1}{dt} &= \frac{1}{\tau} \big( -E_1 + \frac{100(3E_2)_{+}^2}{(120+A_1)^2+(3E_2)_{+}^2} \big) \\ \frac{dE_2}{dt} &= \frac{1}{\tau} \big( -E_2 + \frac{100(3E_1)_{+}^2}{(120+A_1)^2+(3E_1)_{+}^2} \big) \\ \frac{dA_1}{dt} &= \frac{1}{\tau_a} (-A_1 + 0.7 E_1) \\ \frac{dA_2}{dt} &= \frac{1}{\tau_a} (-A_2 + 0.7 E_2) \end{align*}\end{split}\]

where \(\tau = 20\) ms, \(\tau_a = 4000\) ms.

for plotting the isoclines, we consider the first two equations and consider \(A_1\) as a parameter.

def adaptation_isoclines(A1:float, ax:plt.Axes, title:str=""):
    tau = 20
    E1, E2 = sp.symbols('E1 E2')
    X = sp.Matrix([E1, E2])
    F1 = 1/tau *(-E1 + (100 * (3*E2)**2)/ ((120+A1)**2 + (3*E2)**2))
    F2 = 1/tau * (-E2 + (100 * (3*E1)**2)/ ((120+A1)**2 + (3*E1)**2))
    F = sp.Matrix([F1, F2])

    plot_nullclines(F1, F2, symbols=['E1', 'E2'],
                    x_range=(-10,100), y_range=(-10,100),
                    num_points=100, ax=ax, loc='best', title=title)
    eq_points = equilibrium_points(F, X)
    print(eq_points)



fig, ax = plt.subplots(1,2, figsize=(10,4))
adaptation_isoclines(24.0, ax[0], title=r'$A_1=24$')
adaptation_isoclines(36.0, ax[1], title=r'$A_1=36$')

[(0.0, 0.0), (36.0000000000000, 36.0000000000000), (64.0000000000000, 64.0000000000000), (-9.3628088426528 - 18.5411985061749*I, -9.3628088426528 + 18.5411985061749*I), (-9.3628088426528 + 18.5411985061749*I, -9.3628088426528 - 18.5411985061749*I)]
[(0.0, 0.0), (-10.6423173803526 - 21.500641960307*I, -10.6423173803526 + 21.500641960307*I), (-10.6423173803526 + 21.500641960307*I, -10.6423173803526 - 21.500641960307*I), (50.0 - 14.2828568570857*I, 50.0 - 14.2828568570857*I), (50.0 + 14.2828568570857*I, 50.0 + 14.2828568570857*I)]

def ode(x, t, tau, tau_a, stm):
    E1, E2, A1, A2 = x
    dE1dt = 1/tau *(-E1 + (100 * (3*E2+stm*(t>=200)*(t<=400))**2)/ ((120+A1)**2 + (3*E2+stm*(t>=200)*(t<=400))**2))
    dE2dt = 1/tau * (-E2 + (100 * (3*E1)**2)/ ((120+A1)**2 + (3*E1)**2))
    dA1dt = 1/tau_a * (-A1 + 0.7 * E1)
    dA2dt = 1/tau_a * (-A2 + 0.7 * E2)
    return np.array([dE1dt, dE2dt, dA1dt, dA2dt])

t_range = np.arange(0,8000,1)
values = odeint(ode, y0=np.array([1,1,0,0]),
                t=t_range, args=(20, 4000, 50))


fig, ax = plt.subplots(1, figsize=(6,4))
ax.plot(t_range, values[:,0], label="E(t)", c='k', lw=2)
ax.plot(t_range, values[:,2], label='A(t)', c='k', lw=2, ls='--')
ax.legend()
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Spike Rate')
ax.margins(x=0.01)
plt.show()

Competition and neural decisions

\[\begin{split}\begin{align*} \frac{dE_1}{dt} &= \frac{1}{\tau}(-E_1 + S(K_1 - 3E_2)) \\ \frac{dE_2}{dt} &= \frac{1}{\tau}(-E_2 + S(K_2 - 3E_1)) \\ S(x) &= \begin{cases} \frac{100(x)^2}{120^2+(x)^2} & x \geq 0 \\ 0 & x<0 \end{cases} \end{align*}\end{split}\]

def neuron_system(X, t, K1, K2, Tau):
    E1, E2 = X  # Unpack the state vector (E1 and E2 represent the neuron states)

    # Post-synaptic potentials
    PSP1 = (K1 - 3 * E2) * (E2 < K1 / 3)
    PSP2 = (K2 - 3 * E1) * (E1 < K2 / 3)

    # Define the differential equations for each neuron
    dE1_dt = (-E1 + 100 * (PSP1) ** 2 / (120 ** 2 + (PSP1) ** 2)) / Tau
    dE2_dt = (-E2 + 100 * (PSP2) ** 2 / (120 ** 2 + (PSP2) ** 2)) / Tau

    return np.array([dE1_dt, dE2_dt])


def plot_competition(K1, K2, ax, tau=20, y0=[1,0], t_range=np.arange(0,400,5)):
    values = odeint(neuron_system, y0=y0,
                    t=t_range, args=(K1, K2, tau))
    E1, E2 = values[:,0], values[:,1]
    ax[0].plot(t_range, E1, label=r'$E_1(t)$', c='r', alpha=.7)
    ax[0].plot(t_range, E2, label=r'$E_2(t)$', c='b', alpha=.7)
    ax[0].legend()
    ax[0].set_xlabel('Time (ms)')
    ax[0].set_ylabel('Spike Rate')
    ax[0].margins(x=0.01)
    ax[0].set_title(f'$K_1=${K1}, $K_2=${K2}', fontsize=12)

    Xiso = np.arange(0, 101)  # X for Isoclines
    Isocline1 = 100 * (K2 - 3 * Xiso)**2 / (120**2 + (K2 - 3 * Xiso)**2) * (3 * Xiso < K2)
    Isocline2 = 100 * (K1 - 3 * Xiso)**2 / (120**2 + (K1 - 3 * Xiso)**2) * (3 * Xiso < K1)

    ax[1].plot(Xiso, Isocline1, label=r'$dE_1/dt$', c='k', alpha=.7)
    ax[1].plot(Isocline2, Xiso, label=r'$dE_2/dt$', c='k', alpha=.7, ls='--')
    ax[1].plot(E1, E2, label=r'$E_1, E_2$', c='r', alpha=.7)
    ax[1].legend()
    ax[1].set_xlabel('E1')
    ax[1].set_ylabel('E2')
    ax[1].margins(x=0.01)
    return ax

fig, ax = plt.subplots(1, 2, figsize=(10,4))
plot_competition(120, 120, ax)
plt.show()