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