# Lecture 11

## Differential Equations III:

### Steady State Solutions and Mass Action


```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

%matplotlib inline

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

global_fig_width = 10
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
##################################################
```

## Steady-state solutions

- It is often useful to find the **steady state** of a system described by 
  ordinary differential equations


- This occurs when all rates (i.e. derivatives) are zero


- A steady state may be **stable** or **unstable**, depending on whether or not small deviations from the steady state tend to be corrected or amplified


- We can evaluate the stability of a steady state by considering the sign of the derivative nearby

## Example:

Recall the radioactive decay equation $\displaystyle \frac{{\rm d}N}{{\rm d}t} = -\lambda N.$


```python
t = np.linspace(0,6,100)
l = 1.0
N0 = 1.0
def dNdt(N, t):
    return -l * N
N = scipy.integrate.odeint(dNdt, N0, t)

def plot_graph():
    plt.plot(t,N,label='N')
    plt.plot(t,dNdt(N,t),label='dN/dt')
    plt.legend()
    plt.xlabel('t')
```


```python
plot_graph()
```

- Here the steady state occurs when $N=0$, and is 'reached' after infinite time.


- For small positive $N$ the derivative $\;\displaystyle \frac{{\rm d}N}{{\rm d}t}\;$ is negative, moving the 
  system towards the steady state.


- Physical considerations show that negative $N$ is impossible, hence the steady 
  state is stable as we would expect.

## Example:

Consider a simple model of the production and degradation of a protein, shown by the reaction chain $\;\;\overset{v_1}{\rightarrow}S\overset{v_2}{\rightarrow}\;\;$ where $\;v_1\;$ and $\;v_2\;$ are reaction rates, and let $\;s=[S]$.

We can represent the change in concentration of the protein by

$$ \frac{{\rm d}s}{{\rm d}t} = v_1 - v_2 $$

so the steady state occurs when $\;v_1 = v_2$.

### Mass action

$$\;\;\overset{v_1}{\rightarrow}S\overset{v_2}{\rightarrow}\;\;$$

Suppose the reaction is governed by "mass action" kinetics, so $\;v_1=k_1\;$ (constant) and $\;v_2=k_2 s$.

The equation is then

$$\frac{{\rm d}s}{{\rm d}t} = k_1 - k_2 s.$$

The steady state is given by $s = k_1/k_2$, and it is stable:


```python
def plot_linear_steady_state():
    k_1 = 1.0
    k_2 = 1.0
    s = np.linspace(0,2,100)
    dsdt = k_1 - k_2*s
    plt.plot(s,dsdt)
    plt.xlabel('s')
    plt.ylabel(r'$\frac{ds}{dt}$')
    plt.axvline(x=k_1/k_2,ls='--')
    plt.axhline(y=0,ls='--')
```


```python
plot_linear_steady_state()
```

We can solve the differential equation by separation of variables.

$$ \int {1\over k_1 - k_2 s}~{\rm d}s = \int~{\rm d}t $$

$$ s(t) = Be^{-k_2 t} + {k_1\over k_2} $$

Thus the concentration of $S$ relaxes exponentially to the steady state, no matter the initial condition.

If the initial concentration of $S$ is given by $s_0$ then

$$ s(t) = s_0 e^{-k_2 t} + {k_1\over k_2}\left(1-e^{-k_2 t}\right) $$

Suppose instead that $S$ enhances its own production (positive feedback) and the degradation rate is nonlinear, so $v_1=k_1 s$ and $v_2=k_2 s^2$:


$$\frac{{\rm d}s}{{\rm d}t} = k_1s - k_2 s^2$$


Then the system will be at steady state when $k_1 s = k_2 s^2$, i.e. $s=0$ or 
$s=k_1/k_2$.

We can plot $\;\displaystyle \frac{{\rm d}s}{{\rm d}t}\;$ versus $\;s\;$ to see if they are stable or unstable:


```python
def plot_dsdt(k_1,k_2):
    s = np.linspace(-0.1,2,100)
    dsdt = k_1*s - k_2*s**2
    plt.plot(s,dsdt)
    plt.xlabel('s')
    plt.ylabel(r'$\frac{ds}{dt}$')
    plt.axvline(x=k_1/k_2,ls='--')
    plt.axvline(x=0,ls='--')
    plt.axhline(y=0,ls='--')

def plot_quadratic_steady_state():
    interact(plot_dsdt, k_1 = FloatSlider(value=1,min=0.01,max=2.0,step=0.01,continuous_update=False), k_2 = FloatSlider(value=1,min=0.01,max=2.0,step=0.01,continuous_update=False), continuous_update=False)
```


```python
plot_quadratic_steady_state()
```

### Non-graphical method

$$\frac{{\rm d}s}{{\rm d}t} = k_1s - k_2 s^2$$

- When $s$ is small (and positive, since it is a concentration) $s\gg s^2$ so the derivative $k_1 s - k_2 s^2$ will be positive, making $s=0$ an unstable steady state.

- For the other steady state, consider small deviations by $\varepsilon$ from the steady state, e.g. $s=k_1/k_2 + \varepsilon$.
  
Then $\;\;\displaystyle s^2=\frac{k_1^2}{k_2^2} + 2 \varepsilon\frac{k_1}{k_2} + \varepsilon^2\;$ and the derivative is $\; - k_1\varepsilon - k_2\varepsilon^2$.

It is, since $\varepsilon\gg\varepsilon^2$, negative, pushing $s$ back towards the steady state, hence it is a stable steady state.

## Example:

The growth of a cell colony can be modelled by the *logistic* equation

$$ \frac{{\rm d}N}{{\rm d}t} = rN\left(1 - {N\over K}\right) $$

where $\;N(t)\;$ is the number of cells at time $\;t,\;$ and $\;r\;$ and $\;K\;$ are constant parameters (both positive).

The steady state for the system, or equilibrium population size, occurs when the growth rate is zero, i.e.

$$ \frac{{\rm d}N}{{\rm d}t} = 0 \quad\Rightarrow\quad rN = {r\over K}N^2 \quad\Rightarrow\quad N = 0 {\rm~or~} N = K.$$


```python
def plot_dNdt(r,K):
    N = np.linspace(-0.1,2,100)
    dNdt = r*N*(1 - N/K)
    plt.plot(N,dNdt)
    plt.xlabel('N')
    plt.ylabel(r'$\frac{dN}{dt}$')
    plt.axvline(x=K,ls='--')
    plt.axvline(x=0,ls='--')
    plt.axhline(y=0,ls='--')
    
def plot_logistic_steady_state():
    interact(plot_dNdt, r = FloatSlider(value=1,min=0.01,max=2.0,step=0.01),
                    K = FloatSlider(value=1,min=0.01,max=2.0,step=0.01))
```


```python
plot_logistic_steady_state()
```

- for small positive $N$, $rN>0$ and $N<K$ so $\displaystyle \frac{{\rm d}N}{{\rm d}t}>0$ and the population 
  size will increase, meaning that $N=0$ is an unstable steady state.


- In fact the growth rate is positive for $0<N<K$ and negative for $N>K$, making 
  $N=K$ a stable steady state.

We can solve the differential equation to examine the transient behaviour.
This is a separable equation, so

$$ \int {1\over rKN-rN^2}~{\rm d}N = \int {1\over K}~{\rm d}t $$

Solve this, using partial fractions on the left hand side:

- With initial conditions $N(0)=N_0$ you get
$$ N(t) = {K N_0 e^{rt}\over K + N_0(e^{rt}-1)} $$

## Summary of first order differential equations

To solve a differential equation:

1. Calculate the general solution
    1. Try to write it as a separable equation first
    2. Other methods (e.g. integrating factors) not covered in this course
2. This general solution will include an arbitrary constant this may be eliminated using initial conditions (if these are given)
3. Can check your solution numerically using Python

To find steady state (equilibrium) solutions, find points where all (first) derivatives are zero.

To determine stability of these steady states, look at the behaviour of the first derivative in the vicinity of the steady state.
You can sketch the first derivative, or use Python to help with this.

## Chemical reactions to equations: Mass action

- Chemical reactions assume uniform mixing

- Rates proportional to product of masses

- Predator/prey and epidemics have similar rules

1. A and B produce C. The reaction is governed by A meeting B:

$$A+B\underset{^k}{\rightarrow} C$$
    
\begin{align*}
\frac{{\rm d}[A]}{{\rm d}t} &= -k[A][B] \\
\frac{{\rm d}[B]}{{\rm d}t} &= -k[A][B] \\
\frac{{\rm d}[C]}{{\rm d}t} &= k[A][B] \\
\end{align*}

2. Predation of R by W 

$$R+W\underset{^k}{\rightarrow} W$$

$$\frac{{\rm d}R}{{\rm d}t} = -kRW$$

3. Constant production (zeroth order)

$$0 \underset{^k}{\rightarrow} A$$ 

$$\frac{{\rm d}[A]}{{\rm d}t} = k$$

4. Degradation/death 

$$C \underset{^k}{\rightarrow} 0$$

$$\frac{{\rm d}[C]}{{\rm d}t} = -k[C]$$

## Numerically solving differential equations

What if we can't solve the differential equation?

$$ \frac{{\rm d}y}{{\rm d}t} = \sin\left(\frac{y\ln(t\cos(y))}{\sqrt{1+y/e^t}}\right) + 3$$

## Numerically solving differential equations

What if we can't solve the differential equation?

$$ \frac{{\rm d}y}{{\rm d}t} = \sin\left(\frac{y\ln(t\cos(y))}{\sqrt{1+y/e^t}}\right) + 3$$

(or we just don't want to)

## Euler's method

Given a differential equation

$$ \frac{{\rm d}y}{{\rm d}t} = f(y, t)$$

with initial state $\;y(t = t_0) = y_0,\;$ we can *approximate* the state at $t = t_0 + \delta{t}$ as:

$$ y_1 = y(t + \delta{t}) \approx y_0 + f(y, t) \cdot \delta{t}$$

and the next state as
$$ y_2 = y(t + 2\delta{t}) \approx y_1 + f(y, t) \cdot \delta{t}$$

and so on!

## Euler's method

This mean's we can estimate the *entire time course of $y(t)$*, provided: 

1. We can calculate $f(y, t)$ (or approximate it with a computer)
2. We're patient enough to take really tiny steps $\delta{t}$

This is really simple with a computer! See the Python helpsheet for details.

f8f9da05365669fb3633cde46d0aed44961d7ef7
ipynb
Jupyter Notebook
lectures/lecture-11-differential-eq3.ipynb
SABS-R3/2020-essential-maths
5a53d60f1e8fdc04b7bb097ec15800a89f67a047
Apache-2.0
2021-11-27T12:07:13.000Z
2021-11-27T12:07:13.000Z
lectures/lecture-11-differential-eq3.ipynb
SABS-R3/2021-essential-maths
8a81449928e602b51a4a4172afbcd70a02e468b8
Apache-2.0
lectures/lecture-11-differential-eq3.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