<jupyter_start><jupyter_text>Lecture 15 Systems of Differential Equations II: System Simplification<jupyter_code>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
##################################################<jupyter_output><empty_output><jupyter_text>Plan for this lectureRecap- 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 --- RescalingThis 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 forscaling $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 theparameter $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 --- ConservationConsider a chemical system where a product is formed from a substratevia an intermediate complex involving an enzyme.$$S + E\overset{k_1}{\underset{k_{-1}}\rightleftharpoons}C\overset{k_2}\rightarrowP + 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 thenreleased.  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 ODEsystem:\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 stateIn the previous example it may be known that the rate of conversion tocomplex (from $k_1$ and $k_{-1}$) is much faster than the degradationof 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 veryslowly.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 givesan **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 gofrom a 4-dimensional system $(s,e,c,p)$ to a two-dimensional approximationwhich captures some of the behaviour.Two dimensions are good because we can plot their behaviour on a phaseplane diagram...  Phase planes and nullclinesA system of **nonlinear** ODEs may have more than one  fixedpoint (or may have none).  Finding fixed points in two-dimensionalsystems is aided by **nullclines**.A $x$-nullcline is a line where $\dot{x}=0$ and  a $y$-nullcline is aline 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 beannotated to give the direction (and magnitude) of the non-zero derivative.<jupyter_code>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')
plot_nullclines()
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')
plot_phase_plane()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20130010834236187, 'fraction_numerical': 0.031960996749729145, 'mean_word_length': 3.6810344827586206, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '536073', 'license': 'Apache-2.0', 'n_tokens_mistral': 3561, 'n_tokens_neox': 3357, 'n_words': 894, 'path': 'lectures/lecture-15-systems2.ipynb', 'repo_name': 'SABS-R3/2020-essential-maths', 'size': '9230'}