Introduction

What is an ODE

Differential equations can be used to describe the time-dependent behaviour of a variable.
$$\frac{\text{d}\vec{x}}{\text{d}t} = f(\vec{x}, t)$$
The variable stands for a concentration or a number of individuals in a population.

In general, a first order ODE has two parts, the increasing (birth, production,...) and the decreasing (death, degradation, consumption,...) part:

$$\frac{\text{d}\vec{x}}{\text{d}t} = \sum_{}\text{Rates}_{\text{production}} - \sum_{}\text{Rates}_{\text{loss}}$$

You probably already know ways to solve a differential equation algebraically by 'separation of variables' (Trennung der Variablen) in the homogeneous case or 'variation of parameters' (Variation der Konstanten) in the inhomogeneous case. Here, we want to discuss the use of numerical methods to solve your ODE system.

Numerical integration

In principle, every numerical procedure to solve an ODE is based on the so-called "Euler" method. It's very easy to understand, you just have to read the $\frac{\text{d}\vec{x}}{\text{d}t}$ as a $\frac{\Delta \vec{x}}{\Delta t}$. Then you can multiply both sides of the equation with $\Delta t$ and you have an equation describing the change of your variables during a certain time intervall $\Delta t$:

$$ \Delta \vec{x} = f(\vec{x}, t)\cdot \Delta t$$

Of course, the smaller yoy choose the time intervall $\Delta t$, the more accurate your result will be in comparison to the analytical solution.
So it's clear, we chose a tiny one, right? Well, not exactly, the smaller your time intervall the longer the simulation will take. Therefore, we need a compromise and here the provided software will help us by constantly testing and observing the numerical solution and adapt the "step size" $\Delta t$ automatically.

Lotka-Volterra: A prey-predator model

Model equations:

$$ \frac{\mathrm{d}\,R}{\mathrm{d}\,t} = aR-bRW $$$$ \frac{\mathrm{d}\,W}{\mathrm{d}\,t} = cWR-dW $$

Variables:

R: Rabbit population

W: Wolf population

Parameters:

a: rabbit's birth rate

b: predation rate

c: wolf's benefit

d: wolf's death rate

Let's start

we write a small function f, that receives a list of the current values of our variables x, the current time t and parameters p. The function has to evaluate the equations of our system or $\frac{\text{d}\vec{x}}{\text{d}t}$, respectively. Afterwards, it returns the values of the equations as another list.
Important
Since this function f is used by the solver, we are not allowed to change the input (arguments) or output (return value) of this function.

Before we start the simulation of our model, we have to define our system.
We start with our static information:

1. Initial conditions for our variables
2. Values of the paramters
3. Simulation time and number of time points at which we want to have the values for our variables (the time grid). Use numpy!!

Last but not least, we need to import and call our solver. The result will be a matrix with our time courses as columns and the values at the specified time points. Since we have a values for every time point and every species, we can directly plot the results using matplotlib.

Plot results

Phase space diagram

Exercises

1. Implement the Goldbeter model:

$\frac{\mathrm{d}C}{\mathrm{d}t}=v_i-v_dX\frac{C}{K_d+C}-k_dC$

$\frac{\mathrm{d}M}{\mathrm{d}t}=V_1\frac{1-M}{K_1+(1-M)}-V_2\frac{M}{K_2+M}$

$\frac{\mathrm{d}X}{\mathrm{d}t}=V_3\frac{1-X}{K_3+(1-X)}-V_4\frac{X}{K_4+X}$

with $V_1=\frac{C}{K_c+C}V_{M1}$ and $V_3=MV_{M3}$

and $M+M^*=1$ and $X+X^*=1$

Tasks:

• Plot the trajectories and the phase space diagrams
• Perform parameter sampling (20 random parameter sets for the interval [0, 1] ) and plot trajectories and phase space diagrams for each parameter set
• Scan the parameter 'Kd' in a range where the solution oscillates (~20 values) and plot trajectories and phase space diagrams for each parameter set
• Determine frequency and amplitude of those solutions and plot them in dependency of the 'Kd'

Solution