<a href="https://colab.research.google.com/github/cohmathonc/biosci670/blob/master/GrowthModels/GrowthModels_GompertzGrowthNumerical.ipynb" target="_parent"></a>


```python
import numpy as np
import matplotlib.pylab as plt
import pandas as pd
```

# Gompertzian Growth

The Gompertz model is an empirical growth model, based on the observation that population decline can be faster than a decaying exponential, which implies that the death rate must be increasing with age (because a constant death rate would result in exponential decrease).

Mathematically, the Gompertz model can be expressed as a system of two coupled differential equations

\begin{align}
      \frac{dN}{dt} &=G(t)N(t) \tag{1}\\
      \frac{dG}{dt} &=-r G(t) \tag{2}
\end{align}

where
 - the population $N(t)$ grows or shrinks over time with time-dependent growth factor $G(t)$
 - the growth factor $G(t)$ decays exponentially over time with constant rate $r$. 

Here, we explore two ways to solve Gompertz model numerically: 

In the first approach, we make use of the known solution to ODE (2) to solve ODE (1) as a function of population $N$ and time $t$.
Then, we solve the system of coupled ODEs directly, without relying on partial analytic solutions of the problem.

## Use known analytic expression of G(t) and solve single ODE 

We already know the solution to ODE (2):

$$G(t)=G_0\exp(-r\,t)$$

with initial condition $G_0$.

This allows us to rewrite ODE (1) by substituting the functional form of $G(t)$:
$$\frac{dN}{dt}=G(t)N(t) = G_0\exp(-r\,t)\, N(t) \tag{3}$$

Now we have reduced the system of coupled ODEs that define the Gompertz model to a single ODE that depends on time $t$ and population $N$.

This problem can be solved with the tools that we are already familiar with.


```python
# Gompertz ODE with growth rate r and initial growth factor G0
def gompertz_single_ode(t, N):
    r  = 0.3
    G0 = 1
    return G0*np.exp(-r*t)*N

def solve_euler(f, t, y_0):
    y = np.zeros(t.shape)
    y[0] = y_0
    for i in range(0, len(t)-1):
        y[i+1] = y[i] + (t[i+1]-t[i]) * f(t[i],y[i])
    return y

t_single_ode = np.arange(0, 20, 0.1)
y_numeric_single_ode = solve_euler(gompertz_single_ode, t_single_ode, y_0=100)

plt.plot(t_single_ode, y_numeric_single_ode,'b' ,label="numeric solution, single ODE")
plt.legend()
```

## Solve System of Coupled ODEs 

Alternatively, we can solve the system of coupled differential equations numerically:

\begin{align}
      \frac{dN}{dt} &=G(t)N(t) &=& &f(t, N, G) &=& N'\\
      \frac{dG}{dt} &=-r G(t)  &=& &g(t, N, G)     &=& G'
\end{align}

The basic idea of the Euler method was to approximate an unknown function at the next evaluation time point using its current value and its derivative.
Here, we need to approximate two functions simultaneously in each integration step $i$: 

\begin{align}
      N_{i+1} &= N_i + \Delta_i\, f(t_i, N_i, G_i)\\
      G_{i+1} &= G_i + \Delta_i\, g(t_i, N_i, G_i)\\
\end{align}


For implementation, we need a computational function that returns the values $f(t, N, G)$ and $g(t, N, G)$ for given $t, N, G$.
Ideally, we would like this function to use the same call signature as the ODE defining functions that we had used before, i.e. one argument for function values ($N, G$), and one for time.

This can be achieved by vectorization:

Instead of accepting only a single scalar as input for the  current function value, we permit a k-dimensional array as input argument ($k=2$ in this example, i.e. `[N, G]`). 
Likewise, instead of returning only the derivative $y'$, we return an array of derivative values, one for each ODE in our system of ODEs (here: `[N', G']`).


---
**Exercise:**

1. Define a vectorized ODE function with call signature `f(t, y)`, where `t` is an array of evaluation (time) points values and `y` an array of function values. 
 The length of array `y` will correspond to the number of ODEs ($k$) in the coupled system.
 The ODE defining function should return an array `y_p` of same length as `y`.

2. Adapt your implementation of the Euler method to be able to work with this vectorized ODE function. 
   The initial value argument `y_0` will now take an array of $k$ initial values, one for each of the ODEs. The solver function should return an array that contains the solution of all $k$ ODEs at all $n$ integration steps.
   
3. Solve for $N(t)$ and $G(t)$ using initial values $G_0=1$, $N_0=100$ and parameter $r=0.3$ between $t=0$ and $t=20$.
 
4. Plot your results and compare to the result of the previous solution approach, and the analytic solution:

    $$N(t)=N_0e^{-\frac{G_0}{r}(e^{-r\,t} -1)}$$
5. Find a suitable integration stepsize.

6. Verify that the system reaches its carrying capacity $K$ for $t\to+\infty$ and that your solution fulfills the following relation:

    $$N_0e^{\frac{G_0}{r}}=K$$
---


```python
# identify y[0] with N and y[1] with G
def gompertz_coupled_ode(t, y):
    r = 0.3
    y_p = np.zeros(y.shape)
    y_p[0] = y[0] * y[1]
    y_p[1] = -r*y[1]
    return y_p

def solve_euler_ndim(f, t, y_0):
    y = np.zeros((len(y_0),len(t)))
    y[:,0] = y_0
    for i in range(0, len(t)-1):
        y[:,i+1] = y[:,i] + (t[i+1]-t[i]) * f(t[i], y[:,i])
    return y

# initial values
G0 = 1
N0 = 100
# time
t_coupled_ode = np.arange(0, 20, 0.1)

y_numeric_coupled_ode = solve_euler_ndim(gompertz_coupled_ode, t_coupled_ode, y_0=np.array([N0,G0]))
#y_numeric_coupled_ode
```


```python
# plotting
r=0.3
analytic = N0*np.exp(-G0/r * (np.exp(-r*t_coupled_ode)-1))
K = np.exp(G0/r)*N0

fig, axes = plt.subplots(2, 1, figsize=(10,8), sharex=True)

axes[0].plot(t_coupled_ode, y_numeric_coupled_ode[0,:], 'g', label='N(t) - coupled ODE')
axes[0].plot(t_single_ode, y_numeric_single_ode, 'b', label='N(t) - single ODE')
axes[0].plot(t_single_ode, analytic, 'k--', label='N(t) - analytic')
axes[0].axhline(y=K, color='r', linestyle=':', label='expected carrying capacity')
axes[0].set_xlabel("t")
axes[0].set_ylabel("N(t)")
axes[0].legend()

axes[1].plot(t_coupled_ode, y_numeric_coupled_ode[1,:], 'm', label='G(t) - coupled ODE')
axes[1].set_xlabel("t")
axes[1].set_ylabel("G(t)")
axes[1].legend()
plt.show()
```

## Solve Systems of Coupled ODEs with `scipy.integrate.solve_ivp`

`solve_ivp` solves coupled systems *by default*, and always expects an array of initial values `y0` as input.

Applied to the Gompertz model:



```python
from scipy.integrate import solve_ivp

# as in 1.2
def gompertz_coupled_ode(t, y):
    r = 0.3
    y_p = np.zeros(y.shape)
    y_p[0] = y[0] * y[1]
    y_p[1] = -r*y[1]
    return y_p

# initial values
G0 = 1
N0 = 100
# time
t_0 = 0
t_N = 20

sol = solve_ivp(fun=gompertz_coupled_ode, 
                t_span=[t_0, t_N], y0=[N0, G0], 
                rtol=1E-6)   
```


```python
# plotting
r=0.3
analytic = N0*np.exp(-G0/r * (np.exp(-r*t_coupled_ode)-1))
K = np.exp(G0/r)*N0

fig, axes = plt.subplots(2, 1, figsize=(10,8), sharex=True)

axes[0].plot(sol.t, sol.y[0,:], 'g', label='N(t) - coupled ODE (solve_ivp)')
axes[0].plot(t_coupled_ode, analytic, 'k--', label='N(t) - analytic')
axes[0].axhline(y=K, color='r', linestyle=':', label='expected carrying capacity')

axes[0].set_xlabel("t")
axes[0].set_ylabel("N(t)")
axes[0].legend()

axes[1].plot(sol.t, sol.y[1,:], 'm', label='G(t) - coupled ODE')
axes[1].set_xlabel("t")
axes[1].set_ylabel("G(t)")
axes[1].legend()
plt.show()
```

###### About 
This notebook is part of the *biosci670* course on *Mathematical Modeling and Methods for Biomedical Science*.
See https://github.com/cohmathonc/biosci670 for more information and material.

336e2975bb9eb67b2c2f7af32f1b32010a5e1e3c
ipynb
Jupyter Notebook
GrowthModels/GrowthModels_GompertzGrowthNumerical.ipynb
sbranciamore/biosci670coh
92e91aa1e0266b54b0831537b857796eaf0b949c
CC-BY-4.0
2019-02-19T21:20:42.000Z
2021-01-12T20:25:42.000Z
GrowthModels/GrowthModels_GompertzGrowthNumerical.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
GrowthModels/GrowthModels_GompertzGrowthNumerical.ipynb
danielabler/biosci670
a9fc3e3062cc7ee4f2da3c26a4c5fc22c5cb8cb4
CC-BY-4.0
2019-05-15T17:33:41.000Z
2021-02-22T18:11:37.000Z
Qwen/Qwen-72B
1. YES
2. YES