# Lecture 15

## Systems of Differential Equations II:

### System Simplification


```python
import numpy as np
import sympy as sp
import scipy.integrate
sp.init_printing()
##################################################
##### Matplotlib boilerplate for consistency #####
##################################################
from ipywidgets import interact
from ipywidgets import FloatSlider
from matplotlib import pyplot as plt

import cmath

%matplotlib inline

from IPython.display import set_matplotlib_formats
set_matplotlib_formats('svg')

global_fig_width = 8
global_fig_height = global_fig_width / 1.61803399
font_size = 12

plt.rcParams['axes.axisbelow'] = True
plt.rcParams['axes.edgecolor'] = '0.8'
plt.rcParams['axes.grid'] = True
plt.rcParams['axes.labelpad'] = 8
plt.rcParams['axes.linewidth'] = 2
plt.rcParams['axes.titlepad'] = 16.0
plt.rcParams['axes.titlesize'] = font_size * 1.4
plt.rcParams['figure.figsize'] = (global_fig_width, global_fig_height)
plt.rcParams['font.sans-serif'] = ['Computer Modern Sans Serif', 'DejaVu Sans', 'sans-serif']
plt.rcParams['font.size'] = font_size
plt.rcParams['grid.color'] = '0.8'
plt.rcParams['grid.linestyle'] = 'dashed'
plt.rcParams['grid.linewidth'] = 2
plt.rcParams['lines.dash_capstyle'] = 'round'
plt.rcParams['lines.dashed_pattern'] = [1, 4]
plt.rcParams['xtick.labelsize'] = font_size
plt.rcParams['xtick.major.pad'] = 4
plt.rcParams['xtick.major.size'] = 0
plt.rcParams['ytick.labelsize'] = font_size
plt.rcParams['ytick.major.pad'] = 4
plt.rcParams['ytick.major.size'] = 0
##################################################
```

## Plan for this lecture
Recap

- So far we have looked at systems of **first order**, **linear** ODEs in **two dimensions**
    - Are you comfortable with the words in bold?

- These systems can be solved analytically
    - 3 methods of solving systems of first order ODEs
    - Diagonalisation extends to $N$-dimensional systems

Plan

- Aim to look at systems of **first order**, **nonlinear** ODEs in **more dimensions**
- How we go about modelling a problem
- Simplifying systems of ODEs 
    - Reducing number of parameters
    - Reducing number of equations

- Phase plane analysis

## System simplification --- Rescaling


This predator-prey model is in dimensional terms and contains 5 parameters.

\begin{eqnarray*}
\frac{{\rm d}N}{{\rm d}T} &=& rN\left(1-\frac{N}{K}\right) - bNP,\\
\frac{{\rm d}P}{{\rm d}T} &=& ebNP - mP.
\end{eqnarray*}

- $r$: prey intrinsic growth rate
- $K$: prey carrying capacity (bound to logistic growth) 
- $b$: predation rate 
- $e$: conversion efficiency 
- $m$: per capita predator mortality rate

We substitute

- $N=\theta{n}$
- $P=\phi{p}$
- $T=\tau{t}$

where $n,p,t$ are nondimensional, and the $\theta$ is a typical value for $N$ etc.

Note that $\frac{{\rm d}}{{\rm d}T} = \frac{{\rm d}}{{\rm d}t}\frac{{\rm d}t}{{\rm d}T} = \frac{1}{\tau}  \frac{{\rm d}}{{\rm d}t}.$

Original:

\begin{eqnarray*}
\frac{{\rm dN}}{{\rm dT}} &=& rN\left(1-\frac{N}{K}\right) - bNP,\\
\frac{{\rm dP}}{{\rm dT}} &=& ebNP - mP.
\end{eqnarray*}


Rescaled:

\begin{eqnarray*}
\frac{\theta}{\tau}\frac{{\rm d}n}{{\rm d}t} &=& {\theta}rn\left(1-\frac{{\theta}n}{K}\right) - bnp{\theta}{\phi},\\
\frac{\phi}{\tau}\frac{{\rm d}p}{{\rm d}t} &=& ebnp{\theta}{\phi} - mp{\phi}.
\end{eqnarray*}


Note that $\theta,$ $\phi,$ and $\tau$  are arbitrary values for
scaling $N,$ $P,$ and $T$. 


\begin{eqnarray*}
\frac{{\rm d}n}{{\rm d}t} &=& \tau\left( rn\left(1-\frac{{\theta}n}{K}\right) - bnp{\phi}\right),\\
\frac{{\rm d}p}{{\rm d}t} &=& \tau\left( ebnp{\theta} - mp\right).
\end{eqnarray*}


After rearranging we note that choosing $\theta=K$ will remove the
parameter $K$ from the equations.  Choosing $\tau = 1/m$ (or $\tau =
1/r$) will remove a further parameter.


\begin{eqnarray*}
\dot{n} &=& \left( \frac{r}{m}n(1-n) - \frac{b{\phi}}{m}np\right),\\
\dot{p} &=& \left( \frac{ebpK}{m}n - p\right).
\end{eqnarray*}


$\phi=m/b$ gives the two parameter system:


\begin{eqnarray*}
\dot{n} &=& \left( \frac{r}{m}n(1-n) - np\right),\\
\dot{p} &=& \left( \frac{ebpK}{m}n - p\right),
\end{eqnarray*}


where we have shifted the groups of parameters to only two places.

We could make it more obvious that there are only two "tuning dials" by writing it as:


\begin{eqnarray*}
\dot{n} &=&  v_1 n(1-n) - np,\\
\dot{p} &=&  p(v_2 n - 1),
\end{eqnarray*}


## System simplification --- Conservation

Consider a chemical system where a product is formed from a substrate
via an intermediate complex involving an enzyme.
$$
S + E
\overset{k_1}{\underset{k_{-1}}\rightleftharpoons}
C
\overset{k_2}\rightarrow
P + E
$$

This system can be written as a system of 4 ODEs for $s, e, c,
p$ which describe the concentrations of $S, E, C$ and $P$.
\begin{eqnarray*}
\dot{s} &=& -k_1 se + k_{-1}c, \\
\dot{e} &=& -k_1 se + k_{-1}c + k_{2}c,\\
\dot{c} &=& k_1 se - k_{-1}c - k_{2}c, \\ 
\dot{p} &=& k_{2}c.
\end{eqnarray*}

However the enzyme is recycled: it is used in the complex and then
released.  This means that 
$e + c = e_{tot}$ where $e_{tot}$ is constant.

Making the substitution $e =  e_{tot} - c$ to eliminate $e$ we arrive at the 3 ODE
system:


\begin{eqnarray}
\dot{s} &=& -k_1 s(e_{tot} - c) + k_{-1}c, \\
\dot{c} &=& k_1 s(e_{tot} - c) - k_{-1}c - k_{2}c,  \label{eq:cdot}\\
\dot{p} &=& k_{2}c.
\end{eqnarray}


## System simplification --- Quasi-steady state

In the previous example it may be known that the rate of conversion to
complex (from $k_1$ and $k_{-1}$) is much faster than the degradation
of complex into product (rate $k_2$).

If this is the case then we can assume that the concentration of $C$
arrives at some level $c^*$ (a quasi-steady state) and then moves very
slowly.

Solve equation on previous slide for $\dot{c}=0$:

\begin{eqnarray*}
\dot{c} &=& k_1 s(e_{tot} - c) - k_{-1}c - k_{2}c = 0, \\
\implies 
c^*&=&\frac{k_1  e_{tot} s}{k_{-1} + k_2   + k_1 s},
\end{eqnarray*}
which is **not a steady-state**: it moves slowly with $s$.

Substituting this quasi-steady state into the remaining equations gives
an **approximate** system of 2 ODEs:


\begin{eqnarray*}
\dot{s} &=& -\frac{k_1 k_2 e_{tot} s}{ k_{-1} + k_2 + k_1 s}, \\
\dot{p} &=& \frac{k_1 k_2 e_{tot} s}{ k_{-1} + k_2 + k_1 s}.
\end{eqnarray*}


This means that we have used conservation and quasi-steady state to go
from a 4-dimensional system $(s,e,c,p)$ to a two-dimensional approximation
which captures some of the behaviour.

Two dimensions are good because we can plot their behaviour on a phase
plane diagram...

## Phase planes and nullclines
A system of **nonlinear** ODEs may have more than one  fixed
point (or may have none).  Finding fixed points in two-dimensional
systems is aided by **nullclines**.

A $x$-nullcline is a line where $\dot{x}=0$ and  a $y$-nullcline is a
line where $\dot{y}=0$. (These are straight lines in the linear case.)

For example


\begin{eqnarray*}
\dot{x} &=& x(1-x) -xy,\\
\dot{y} &=& 2y\left(1-\frac{y}{2}\right) - 3xy,
\end{eqnarray*}


has $x$-nullclines at $x=0$ and $1-x-y=0$; and $y$-nullclines at $y=0$ and
$2-3x-y=0$.

Nullcline intersections give us the fixed points.  Nullclines can be
annotated to give the direction (and magnitude) of the non-zero derivative.


```python
def plot_nullclines():
    
    plt.figure(figsize=(10,10))
    
    # plot nullclines
    x = np.linspace(-0.1,1.1,24)
    y = np.linspace(-0.1,2.1,24)
    
    plt.hlines(0,-1,15, color='#F39200', lw=4, label='y-nullcline 1')
    plt.plot(x,1 - x, color='#0072bd', lw=4, label='x-nullcline 2')
    plt.vlines(0,-1,15, color='#0072bd', lw=4, label='x-nullcline 1')
    plt.plot(x,2 - 3*x, color='#F39200', lw=4, label='y-nullcline 2')

    plt.xlim(-0.05,1.1)
    plt.ylim(-0.05,2.1)
    plt.xlabel('x')
    plt.ylabel('y')
```


```python
plot_nullclines()
```


```python
def dX_dt(X, t):
    return np.array([ X[0]*(1. - X[0]) - X[0]*X[1],
                     2.*X[1]*(1.-X[1]/2.) -3*X[0]*X[1]])

def plot_phase_plane():
    
    plt.figure(figsize=(10,10))
    
    init_x = [1.05, 0.9, 0.7, 0.5, 0.5, 0.32, 0.1]
    init_y = [1.0, 1.3, 1.6, 1.8, 0.2, 0.2, 0.2]
    
    plt.plot(init_x, init_y, 'g*', markersize=20)
    
    for v in zip(init_x,init_y):
        X0 = v                              # starting point
        X = scipy.integrate.odeint( dX_dt, X0, np.linspace(0,10,100))         # we don't need infodict here
        plt.plot( X[:,0], X[:,1], lw=3, color='green')
    
    
    
    # plot nullclines
    x = np.linspace(-0.1,1.1,24)
    y = np.linspace(-0.1,2.1,24)
    
    plt.hlines(0,-1,15, color='#F39200', lw=4, label='y-nullcline 1')
    plt.plot(x,1 - x, color='#0072bd', lw=4, label='x-nullcline 2')
    plt.vlines(0,-1,15, color='#0072bd', lw=4, label='x-nullcline 1')
    plt.plot(x,2 - 3*x, color='#F39200', lw=4, label='y-nullcline 2')

    # quiverplot - define a grid and compute direction at each point
    X , Y  = np.meshgrid(x, y)                  # create a grid
    DX = X*(1-X) - X*Y                          # evaluate dx/dt
    DY = 2*Y*(1 - Y/2.0) - 3*X*Y                # evaluate dy/dt               
    M = (np.hypot(DX, DY))                      # norm growth rate 
    M[ M == 0] = 1.                             # avoid zero division errors 

    plt.quiver(X, Y, DX/M, DY/M, M)
    plt.xlim(-0.05,1.1)
    plt.ylim(-0.05,2.1)
    plt.xlabel('x')
    plt.ylabel('y')
```


```python
plot_phase_plane()
```

## Summary

- Simplification
    - Rescaling to dimensionless quantities
    - Conservation
    - Quasi-steady state approximation
  
- Nullclines for finding steady states and phase flow

5566388a988b2479b795fa4768d9c2f735ee516e
ipynb
Jupyter Notebook
lectures/lecture-15-systems2.ipynb
SABS-R3/2020-essential-maths
5a53d60f1e8fdc04b7bb097ec15800a89f67a047
Apache-2.0
2021-11-27T12:07:13.000Z
2021-11-27T12:07:13.000Z
lectures/lecture-15-systems2.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
lectures/lecture-15-systems2.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
2020-10-30T17:34:52.000Z
2020-10-30T17:34:52.000Z
Qwen/Qwen-72B
1. YES
2. YES