import numpy as np
import random
import matplotlib.pyplot as plt

A dummy example of how to numerically solve a system of first-order ODEs with Euler's method

$$ \begin{align} &\frac{dx}{dt}=+p_1 x + p_2 y - p_3 z\\ &\frac{dy}{dt}=+p_1 x - p_2 y + p_3 z\\ &\frac{dz}{dt}=-p_1 x + p_2 y + p_3 z \end{align} $$

def ode_example(sol0,p,tspan=[0,1],dt=0.01):
    
    # sol0 is the initial state of the system
    
    # vector p contains all the fixed parameters
    
    [p1,p2,p3] = p
    
    # define functions to calculate rate of change for each variable
    
    def dx(x,y,z):
        return p1*x  +p2*y -p3*z
    def dy(x,y,z):
        return p1*x  -p2*y +p3*z
    def dz(x,y,z):
        return -p1*x +p2*y +p3*z
    
    # time points where solution should be calculated
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    
    # initialize the solution arrays
    x = [sol0[0]]
    y = [sol0[1]]
    z = [sol0[2]]
    
    for t in T[0:-1]:
        x.append(x[-1] + dt * dx(x[-1],y[-1],z[-1]))
        y.append(y[-1] + dt * dy(x[-1],y[-1],z[-1]))
        z.append(z[-1] + dt * dz(x[-1],y[-1],z[-1]))
    
    return T,x,y,z

sol0 = [0.3,0.2,0.1] # initial state
tspan = [0,2] # the range of time
dt = 0.001 # time step
p = [-0.1,0.2,0.1] # parameters

T,x,y,z = ode_example(sol0,p,tspan,dt) # solve!

plt.plot(T,x,label="x")
plt.plot(T,y,label="y")
plt.plot(T,z,label="z")
plt.legend()

# plt.plot(u,p) # to plot the trajectory in phase diagram

plt.show()

The PhoP/PhoQ system

$$ \begin{align} &\frac{dU}{dt} = k_1 + k_2\frac{P^n}{K^n+P^n}+k_a^-P-(k_d+k_a^+)U\\ &\frac{dP}{dt}=k_a^+U-(k_d+k_a^-)P \end{align} $$

## function to solve deterministically using Euler's method

def ode(sol0,k,n=1,tspan=[0,1],dt=0.01):
    
    [k1,k2,ka1,ka2,kd,K] = k
    
    def du(u,p):
        return k1 + k2*p**n/(K**n+p**n)+ka1*p-(kd+ka2)*u
    def dp(u,p):
        return ka2*u-(kd+ka1)*p
    
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    u = [sol0[0]]
    p = [sol0[1]]
    
    for t in T[0:-1]:
        u.append(u[-1] + dt * du(u[-1],p[-1]))
        p.append(p[-1] + dt * dp(u[-1],p[-1]))
    
    return T,u,p
    

k1 = 0.01
k2 = 0.3
ka1 =  5
ka2 =  1
kd = 0.08
K = 0.2
n = 3.0
k = [k1,k2,ka1,ka2,kd,K]

## solving the deterministic ODE
sol0 = [0.1,0.1] # lower equilibrium
sol1 = [0.5,0.5] # higher equilibrium
tspan = [0,30]
dt = 0.001
T,u,p = ode(sol1,k,n,tspan,dt)

plt.plot(T,u,label="U")
plt.plot(T,p,label="P")
plt.legend()

# plt.plot(u,p) # to plot the trajectory in phase diagram

plt.show()

## showing the phase diagram

def phase(k,n,xlims=[0,1],ylims=[0,1],dx=0.1,dy=0.1):
    
    [k1,k2,ka1,ka2,kd,K] = k
    
    def du(u,p):
        return k1 + k2*p**n/(K**n+p**n)+ka1*p-(kd+ka2)*u
    def dp(u,p):
        return ka2*u-(kd+ka1)*p
    
    x = np.arange(xlims[0],xlims[-1]+dx,dx)
    y = np.arange(ylims[0],ylims[-1]+dy,dy)
    Y,X = np.mgrid[ylims[0]:ylims[-1]:dy,xlims[0]:xlims[-1]:dx]
    U = du(X,Y)
    P = dp(X,Y)
    
    speed = np.sqrt(U**2 + P**2)
    lw = 5*speed / speed.max()
#     plt.quiver(X,Y,U,P) # alternative method of plotting
    plt.streamplot(X,Y,U,P,density=1,linewidth=lw)
    plt.xlabel("U")
    plt.ylabel("P")

# lower equilibrium

phase(k,n,xlims=[0,0.2],ylims=[0,0.05],dx=0.01,dy=0.01)

plt.show()

# higher equilibrium

phase(k,n,xlims=[0,6],ylims=[0,1.5],dx=0.01,dy=0.01)

plt.show()

Phase diagram of an oscillating system

# https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

# Lotka-Volterra equations

def phase_LV(p,xlims=[0,1],ylims=[0,1],dx=0.1,dy=0.1):
    
    [α,β,γ,δ] = p
    
    def dPrey(x,y):
        return α*x - β*x*y
    def dPredator(x,y):
        return γ*x*y - δ*y

    
    cor_x = np.arange(xlims[0],xlims[-1]+dx,dx)
    cor_y = np.arange(ylims[0],ylims[-1]+dy,dy)
    mesh_y,mesh_x = np.mgrid[ylims[0]:ylims[-1]:dy,xlims[0]:xlims[-1]:dx]
    
    Prey = dPrey(mesh_x,mesh_y)
    Predator = dPredator(mesh_x,mesh_y)
    
    speed = np.sqrt(Prey**2 + Predator**2)
    lw = 5*speed / speed.max()

    plt.streamplot(mesh_x,mesh_y,Prey,Predator,density=1,linewidth=lw)
    plt.xlabel("Prey")
    plt.ylabel("Predator")

α = 1.1
β = 0.4
γ = 0.5
δ = 0.2
p = [α,β,γ,δ]

phase_LV(p,xlims=[0,4],ylims=[0,8],dx=0.01,dy=0.01)

plt.show()

## function to solve deterministically using Euler's method

def ode_LV(sol0,p,tspan=[0,1],dt=0.01):
        
    [α,β,γ,δ] = p
    
    def dPrey(x,y):
        return α*x - β*x*y
    def dPredator(x,y):
        return γ*x*y - δ*y
    
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    prey = [sol0[0]]
    predator = [sol0[1]]
    
    for t in T[0:-1]:
        prey.append(prey[-1] + dt * dPrey(prey[-1],predator[-1]))
        predator.append(predator[-1] + dt * dPredator(prey[-1],predator[-1]))
    
    return T,prey,predator

α = 1.1
β = 0.4
γ = 0.5
δ = 0.2
p = [α,β,γ,δ]

sol0 = [10,10] # one cycle
sol1 = [1,4] # another cycle

tspan = [0,150]
dt = 0.001

T,prey,predator = ode_LV(sol1,p,tspan,dt)

# plt.plot(T,prey,label="Prey")
# plt.plot(T,predator,label="Predator")
# plt.legend()

plt.plot(prey,predator) # to plot the trajectory in phase diagram
plt.xlabel("Prey")
plt.ylabel("Predator")

plt.show()

See animation: https://www.geogebra.org/m/y746ry8g

Stochasticity produces nontrivial behaviors in some systems

In the L-V oscillation model, stochastic deviation from the idealized cycle can have lasting consequences..

The trajectory no longer stays close to a deterministic solution, because all cycles are equally stable. There is no "force" pulling the system back to its deterministic trajectory, so that stochasticity moves the system around different cycles randomly.

def ode_LV_noisy(sol0,p,tspan=[0,1],dt=0.01,σ=0.1):
        
    [α,β,γ,δ] = p
    
    def dPrey(x,y):
        return α*x - β*x*y
    def dPredator(x,y):
        return γ*x*y - δ*y
    
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    prey = [sol0[0]]
    predator = [sol0[1]]
    
    for t in T[0:-1]:
        prey.append(prey[-1] + dt * dPrey(prey[-1],predator[-1]) + np.random.normal(scale=dt * (σ**2) * prey[-1]))
        predator.append(predator[-1] + dt * dPredator(prey[-1],predator[-1]) + np.random.normal(scale=dt * (σ**2) * predator[-1]))
    
    return T,prey,predator

α = 1.1
β = 0.4
γ = 0.5
δ = 0.2
p = [α,β,γ,δ]

sol0 = [10,10] # one cycle
sol1 = [1,4] # another cycle

tspan = [0,1500]
dt = 0.1

T,prey,predator = ode_LV_noisy(sol1,p,tspan,dt,σ=0.01)

# plt.plot(T,prey,label="Prey")
# plt.plot(T,predator,label="Predator")
# plt.legend()

plt.plot(prey,predator) # to plot the trajectory in phase diagram
plt.xlabel("Prey")
plt.ylabel("Predator")

plt.show()

Noise interacts with bistable systems and causes the system to jump between equilibria

1D approximation of U (not covered in section, not required to know!)

If $k_a^+,k_a^-$ are much larger than other parameters, then the conversion between $U$ and $P$ happens at a much shorter time scale than the other reactions. We may then assume that, at the time scale of conversion, all other reactions are too slow to affect the system in any way. This conversion process brings the system to the line of conversion equilibriums $P = k_a^+ U/(k_a^-+k_a^+)$. You can see this line in the previous phase diagram! It is the thin line that starts at 0 and absorbs many thick trajectories.

Once you are on $P = k_a^+ U/(k_a^-+k_a^+)$, the conversion has reached its equilibrium, and it is only then that other slow reactions begin to show their effects on the system. The system will evolve along the line $P = k_a^+ U/(k_a^-+k_a^+)$, governed by the slow reactions.

That's the reason why you see the deterministic solution has two phases. In the first phase, U and P rapidly increase/decrease to a turning point, and after that, both of them change slowly through time.

It is possible to simplify the system further along that line $P = k_a^+ U/(k_a^-+k_a^+)$, and you only have to work with one variables now because $U$ and $P$ are locked together by that relationship. It can be shown that the system evolves following:

$$ \frac{dU}{dt} = \frac{1}{\bar{k}_a+1}\left(k_1+\frac{k_2\bar{k}_a^n U^n}{K^n+\bar{k}_a^nU^n}\right)-k_d U $$

where $\bar{k}_a = k_a^+/k_a^-$

## showing the one-dimensional approximation along the conversion equilibrium line

u = np.arange(0,6,0.01)
ut = (k1 + k2*((ka2/ka1)**n)*(u**n)/(K**n + ((ka2/ka1)**n)*(u**n)))/(1+ka2/ka1) - kd * u
plt.plot(u,ut)
plt.plot(u,np.zeros(len(u)))
plt.title("1D approximation of dU/dt")
plt.xlabel("U")
plt.ylabel("dU/dt")

Text(0, 0.5, 'dU/dt')

def ode2(sol0,k,n=1,tspan=[0,1],dt=0.01):
    
    [k1,k2,ka1,ka2,kd,K] = k
    
    def du(u):
        return (k1 + k2*((ka2/ka1)**n)*(u**n)/(K**n + ((ka2/ka1)**n)*(u**n)))/(ka2/ka1+1) - kd * u
    
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    u = [sol0[0]]
    
    for t in T[0:-1]:
        u.append(u[-1] + dt * du(u[-1]))
    
    return T,u

Here, you can see that the 1-D approximation of U (dotted line) is very accurate, it is only slightly different from the real trajectory (blue line)

## solving the deterministic ODE along with the 1-d approximation

k1 = 0.01
k2 = 0.3
ka1 =  5
ka2 =  1
kd = 0.08
K = 0.2
n = 3.0
k = [k1,k2,ka1,ka2,kd,K]

sol0 = [0.5,0.1]
tspan = [0,150]
dt = 0.001
T,u,p = ode(sol0,k,n,tspan,dt)
T,u2 = ode2([(sol0[0]+sol0[1])*ka1/(ka2+ka1)],k,n,tspan,dt)

plt.plot(T,u,label="U")
plt.plot(T,p,label="P")
plt.plot(T,u2,label="U,approx",linestyle="dotted")
plt.legend()

plt.show()

Next, we add noise to the 1-D equation. Noise is assumed to be Normal with mean 0 and standard deviation $\sqrt{k_1 u \times dt}$

def sde(sol0,k,n=1,tspan=[0,1],dt=0.01):
    
    [k1,k2,ka1,ka2,kd,K] = k
    
    def du(u):
        return (k1 + k2*((ka2/ka1)**n)*(u**n)/(K**n + ((ka2/ka1)**n)*(u**n)))/(ka2/ka1+1) - kd * u
    
    T = np.arange(tspan[0],tspan[-1]+dt,dt)
    u = [sol0[0]]
    
    for t in T[0:-1]:
        u.append(u[-1] + dt * du(u[-1]) + np.random.normal(loc=0.0, scale=np.sqrt(dt)) * np.sqrt(k1 * np.abs(u[-1])) * (u[-1]>0))
    
    return T,u

This figure shows that sometimes noisy system stays near the deterministic equilibrium:

## solving the deterministic ODE along with the 1-d approximation
# k1 = 0.05
# n = 9.0
# ka1 =  20
# ka2 =  5
sol0 = [0.2,0.2]
tspan = [0,1000]
dt = 0.01
# T,u,p = ode(sol0,k,n,tspan,dt)
T,u2 = ode2([(sol0[0]+sol0[1])*ka1/(ka2+ka1)],k,n,tspan,dt)
T,u3 = sde([(sol0[0]+sol0[1])*ka1/(ka2+ka1)],k,n,tspan,dt)

plt.plot(T,u2,label="U,approx",linestyle="dotted")
plt.plot(T,u3,label="U,with fluctuation",linestyle="dotted")
plt.legend()

plt.show()

The following figure shows that noise can push the system to a different equilibrium after some long waiting time. The parameters are exactly the same as the previous one

## solving the deterministic ODE along with the 1-d approximation
# k1 = 0.05
# n = 9.0
# ka1 =  20
# ka2 =  5
sol0 = [0.2,0.2]
tspan = [0,1000]
dt = 0.01
# T,u,p = ode(sol0,k,n,tspan,dt)
T,u2 = ode2([(sol0[0]+sol0[1])*ka1/(ka2+ka1)],k,n,tspan,dt)
T,u3 = sde([(sol0[0]+sol0[1])*ka1/(ka2+ka1)],k,n,tspan,dt)

plt.plot(T,u2,label="U,approx",linestyle="dotted")
plt.plot(T,u3,label="U,with fluctuation",linestyle="dotted")
plt.legend()

plt.show()