4 Differential equations

4.1 Introduction

A modern drug to combat malaria is artemisinine. Artemisinine is degraded in the body. One hour after intake, half the quantity of the drug has disappeared, and an hour later again half of what was remaining, etc.. To estimate the effectiveness and correct dosage of of the drug, and also to be able to synthesize more stable forms of the drug, it is important to be able to understand the dynamics of drug degradation, or its pharmaco-kinetics. The relation between time and artemisinine concentration in the body can be described relatively easily. Namely, at any point in time, the rate at which artemisinine disappears from the body is proportional to the concentration of artemisinine. In mathematical terms: if $q(t)$ is the function that gives the concentration of artemisinine in the body at any point in time $t$, then the change of $q(t)$ in time, $q'(t)$, that is the first order derivative of $q$ to time, equals

\[ q'(t) = - k \, q(t) \tag{4.1}\]

where $k$ is a constant. This relation is called “first-order kinetics” because the change of $q$ is proportional to the first order of $q$ ($q^1=q$). Also, note the minus sign in front of $k$. $k$ Is defined as a positive number, and the minus sign indicates that we are looking at a decrease of $q$.

The equation Equation 4.1 is an example of a differential equation: an equation for the derivative of a function. Differential equations play a key role in biology and other natural sciences, and from the example it is clear why that is the case: when it is possible to describe the rate of change of a variable, but not the variable itself, which is very often the case, you will be dealing with a differential equation.

In this chapter we will discuss some mathematical aspects of differential equations: different types, whether they have a solution, how you find the solution, and —because solving differential equations can be rather cumbersome— what you can still do if you can’t find a solution.

4.2 Differential equations and their solutions

Differential equations differ from “normal” equations (the correct term is algebraic equations) because they contain a derivative; but also because their solutions are not one or a few numbers but complete functions. An algebraic equation like $2x = 4$ is an assertion about a number, which can be true or not; $x = 2$ is called a solution of an equation because the assertion is true for this number: after substitution of $x = 2$ in $2x = 4$ both sides of the equation become equal. In the same manner, a differential equation is an assertion about a function, and a function is the solution of the differential equation when the assertion is true for that function, i.e. both sides of the differential equation become equal after substitution of the solution. As an example we look again at the artemisinine concentration.

Example 4.1 \[ q'(t) = - 0.9 \, q(t) \tag{4.2}\]

We just try a function that decreases with time, for example a parabola in time. The number $-0.9$ will probably also have to be in it, so we use that as a coefficient:

\[ q_{1}(t) = 1 - 0.9 \, t^2 \] Is this a solution for the differential equation Equation 4.2? To answer that question we substitute this trial function for $q(t)$ on both sides of the equation¹. For the left hand side we need the first order derivative $q_{1}'(t)$ of our trial function, which equals $-1.8 \, t$. For the right hand side we get $- 0.9 \, (1 - 0.9 \, t^2) = -0.9 + {(0.9)}^2 \, t^2$. Right and left hand sides are not equal, so our trial function can not be a solution to the differential equation Equation 4.2.

We now try an exponentially decreasing function:

\[ q_{2}(t) = e^{-0.9 \, t} \]

The derivative (left hand side of the differential equation) equals $q_{2}'(t) = -0.9 \, e^{-0.9 \, t}$, and the right hand side of the differential equation becomes $-0.9 e^{-0.9 \, t}$. Left and right hand sides are equal with this function, so the function $q_2$ is a solution to the differential equation! It is not the only solution, however. Also the function

\[ q_{3}(t) = C e^{-0.9 \, t} \] is a solution for any constant $C$ (check that). If we interpret $C$ in terms of the pharmaco-kinetics of artemisinine then this constant equals the concentration of artemisinine at time $t=0$, since $C e^0 = C$. It is not really surprising that this constant can not be determined by the differential equation, since the differential equation only gives us information about the change of the aremisinine concentration, but not about the initial concentration. This is a general property of differential equations: the functions that are solutions to differential equations always contain one or more undetermined constants. In the formal sense, a differential equation therefore always has a set of solutions, in this case $C e^{-0.9 \, t}$, in which all possible values of $C$ are valid.

The example gives rise to a number of questions.

For this differential equation we found a large number of —similar— solutions. Was that all, or are there more?
Was this a special case or does it happen more often that a differential equation has many solutions?
We found the solution by trial. Is there a more systematic way?
How can you know whether a differential equation has any solution at all? If it does not exist you would be looking for one forever!

Fist questions 1 and 2: the many solutions. We saw in the example that every function of the kind $C e^{-0.9 \, t}$ is a solution, but that’s all! The same thing happens with all differential equations. There is never exactly one solution, but a general solution, in which there is an arbitrary constant. If you want to determine this constant, you have to have additional information, apart from the differential equation, about the solution. This extra bit of information is called a boundary condition. In many situations that we will discuss this is often the value of the function at $t=0$, the initial value.

All of this can be expressed exactly in the following proposition. It contains the term order. A differential equation in which only the first derivative $f'$ of a function $f$ appears is a first-order differential equation. If $f''$ also appears in it, it would be called a second-order differential equation, etc.

Theorem 4.1 A differential equation of the $n$-th order to which $n$ boundary conditions are applied has exactly one solution.

Is this proposition valid for all differential equation? Strictly, no. But in practice: yes. The differential equation has to satisfy certain (rather technical) conditions, but these encompass most cases. The cases that you will encounter are very likely to comply. With this, questions 1 and 2 are answered and 4 also.

Question 3 remains open. Because, although theory tells us that a differential equation has a solution, it does not mean that you will find it.

The remainder of this section is concerned with this: when can you find a solution and what can you do if you don’t succeed. In fact, solving a differential equation is not the only thing you can do. To handle differential equations there are roughly three methods:

Solve it
Analyze the equation itself
Computer simulation

We will discuss these methods in the next three paragraphs.

4.3 Classification, handling and solving

We already indicated that solving differential equations may be cumbersome. Although, not always: for some types of differential equations there exist standard methods that yield a solution guaranteed. For others there are methods that sometimes work, sometimes not. Then there is also a category that needs some juggling. For many differential equations it’s just hopeless. As a user, you should have an idea whether a differential equation is easily solvable or not. For that reason we first discuss the classification of differential equations. There is a number of criteria for the classification of differential equations, of which we mention the most important. They are independent, and all combinations occur.

The order of a differential equation

The order of a differential equation can be 1, 2, 3, etc. and is related to the highest derivative appearing in the equation. In biology you will most often find first-order equations, or sometimes second-order.

Partial differential equations

Partial differential equations are distinguished from ordinary differential equations by the fact that they are concerned with functions of two or more variables (often: place and time), and indicate a relation between the (partial) derivatives to these variables. The occur often in biology and are important, but difficult to handle. They will, therefore, not be discussed here.

Linear differential equations

The last category we will discuss are linear differential equations. A differential equation is called linear when the right side is a linear function of the unknown function(s). To illustrate, two differential equations

\[ y'(t) = 3y(t) + \sin{t} \quad \text{(A)} \qquad y'(t) = 2y(t) - 3 y^2(t) \quad \text{(B)} \]

Linearity of a differential equation is determined by the way in which $y(t)$ in the right hand side. Example A is linear, example B is non-linear. The criterion is whether the right hand side of the differential equation can be written as something times $y(t)$ plus something else, which can be done in case A (the term $\sin{t}$ does not matter in this case) but not with B. The only thing that matters is the way in which $y(t)$ occurs in the right hand side.

In addition to individual differential equations, there also exist systems of (coupled) differential equations. These concern two or more functions that determine each other’s derivatives.

Example 4.2 To illustrate, some examples of different types of differential equations.

$\displaystyle y'(t)=ay(t)$	ordinary, first order, linear
$\displaystyle y'(t) = a {\left\{ y(t) \right\}}^{\frac{2}{3}} - b y(t)$	ordinary, first order, nonlinear
$\displaystyle y'(t) = ay(t) + b$	ordinary, first order, linear
$\displaystyle y'(t)=\frac{a y(t)}{b + y(t)}$	ordinary, first order, nonlinear
$\displaystyle \left\{ \begin{array}{l} y'(t)=ay(t) + bz(t)\\ z'(t)=cy(t)-dz(t)\end{array} \right.$	ordinary, first order system, linear
$\displaystyle y'(t) = ay(t) + b\sin{at}$	ordinary, first order, linear (!)
$\displaystyle y''(t) = ay'(t) + by(t)$	ordinary, second order, linear
$\displaystyle \pp{C}{t} = D \left( \frac{\partial^2 C}{\partial x^2} + \frac{\partial^2 C}{\partial y^2} \right)$	partial, second order, linear

We will now discuss a few solutions of common differential equations.

Linear ordinary differential equations

First, linear ordinary differential equations. The simplest and most orderly type, with a nice property: they can always be solved. This makes them very attractive for the modeler. For nonlinear ordinary differential equations the case is often more complicated. The general form of a linear (first order) ordinary differential equation is

\[ y'(t) = f(t)y(t) + g(t) \]

in which both coefficients $f(t)$ and $g(t)$ are allowed to be functions of $t$, but could also be constant. The solution is in principle always an $e$-power. In case that the coefficient are constants the solution is very simple. In the following we show a few solutions from the most simple to the most general case. The boundary condition is $y(0)=y_0$ for all these cases.

\[ \begin{align*} y'(t) & = k y(t) & y(t) & = y_0 e^{kt} \\ y'(t) & = k y(t) + a & y(t) & = \left( \frac{a}{k} + y_0 \right) e^{kt} - \frac{a}{k} \\ y'(t) & = k y(t) + g(t) & y(t) & = e^{kt}\left( \int\limits_0^t g(u)e^{-ku} \; du + y_0 \right) \\ y'(t) & = f(t) y(t) + g(t) & y(t) & = e^{F(t)} \left( \int\limits_0^t g(u)e^{-F(u)} \; du + y_0 \right) \end{align*} \]

\[ \qquad \qquad \text{in which } F(t) = \int\limits_0^t f(v) \; dv \] Although the last two solutions contain integrals, they are solutions: they are integrals of known functions. However, we will not encounter these integrals in this course, and they are only mentioned here for completeness.

Example 4.3 (Pharmacokinetics) Our example \[ q'(t) = -k q(t) \] clearly is a linear equation of the first type in the list above. We concluded that this equation described the degradation of a compound, for example a drug, in the body. The solution of this differential equation is, therefore, \[ q(t) = q_0 e^{-k t} \]

A more interesting equation is one that describes what happens if a drug is administered and degraded at the same time.

Example 4.4 (A drug that is degraded while administered) Suppose that the administration via an infusion is a constant amount $a$ per time unit. The corresponding differential equation will then be

\[ q'(t) = a - k q(t) \tag{4.3}\]

This is also a linear differential equation, in which $g(t)=a$ is a constant function. The solution is given by the second row in the list above:

\[ q(t) = \left( \frac{a}{-k} + q_0 \right) e^{-kt} - \frac{a}{-k} = \left( q_0 - \frac{a}{k} \right) e^{-kt} + \frac{a}{k} \] The course of this function is shown in Figure 4.1}. From this figure, as well as from its equation, you can derive that the concentration of the drug will stabilize at a value $\frac{a}{k}$. A mathematician would say that

\[ \lim_{t \to \infty} \left( q_0 - \frac{a}{k} \right) e^{-kt} + \frac{a}{k} = \frac{a}{k} \]

Figure 4.1: Time course of the concentration, $q(t)$, of a drug that is administered at constant rate $a$ and is degraded with a rate constant $k$ according to Equation 4.3. The boundary condition is $q_0=0$.

Logistic equations

In this section we discuss a common nonlinear ordinary differential equation, the so-called logistic equation, with boundary condition $y(0)=y_0$.

\[ y'(t)=ay(t) \left\{ 1 - \frac{y(t)}{b} \right\} \qquad y(t) = \frac{b}{1 - Ce^{-at}} \quad \text{with} \quad C = \frac{b}{y_0} -1 \tag{4.4}\]

This is the form in which you will encounter it in applications in ecology. The more `neutral’ form is

\[ y'(t) = a y(t) + \tfrac{1}{b} y^2(t) = a y(t) + b' y^2(t) \]

which is fully equivalent (with $b'=\tfrac{1}{b}$).

*Solving a differential equation by separation of variables

It is sometimes possible to solve a differential equation by so-called separation of variables. When a differential equations

\[ \dd{x}{t} = f(x,t) \]

can be written as

\[ \dd{x}{t} = g(x) \cdot h(t) \]

then the solution can be obtained from the equation

\[ \int \frac{1}{g(x)} \; dx = \int h(t) \; dt + C \] if you can rework this to an equation that expresses $x$ explicitly as a function of $t$.

Example 4.5 (Exponential growth or decrease) An example is the well-known differential equation for proportional growth (positive $k$) or decrease (negative $k$) of a variable $x$:

\[ \dd{x}{t} = kx \] of which the right hand side can be written as $g(x) \cdot h(t)$ with $g(x) = x$ and $h(t) = k$. The solution can be obtained from

\[ \int \frac{1}{x} \; dx = \int k \; dt + C \] or

\[ \ln{(x)} = kt + C \] Reworking this equation, we can obtain $x$ as an explicit function of $t$:

\[ x(t) = C' \cdot e^{k t} \] where the constant $C' = e^{C}$.

The rule for the separation of variables can be derived from the chain rule for differentiation of functions. Namely, if the right hand side of the differential equation can be written as $g(x) h(t)$ then

\[ \frac{1}{g(x)}\dd{x}{t} = h(t) \] The left hand side of this equation can be thought of as the derivative of a composite function $u(x(t))$:

\[ \dd{u(x(t))}{t} = \dd{u}{x} \dd{x}{t} \] and the solution we are looking for can be written as $u^{-1}(H(t))$ where $u^{-1}(x)$ is the inverse of the function $u(x)$ which itself is the anti-derivative of $\frac{1}{g(x)}$ and $H(t)$ is the anti-derivative of $h(t)$. This works, provided that we can find these anti-derivatives and that the inverse function of $u(x)$ exists.

Systems of linear ordinary differential equations

Finally, we will discuss systems of two linear ordinary differential equations, with constant coefficients. Each of both solutions is a sum of two $e$-powers. First we discuss a very common special case, and then the more general and complicated case.

Boundary conditions: $x(0) = x_0$, $y(0) = y_0$.

\[ \left\{ \begin{array}{l} x'(t)=-ax(t) \\ y'(t)=bx(t)-cy(t) \end{array} \right. \; \text{sol:} \; \left\{ \begin{array}{l} \displaystyle x(t) = x_0e^{-at} \\ \displaystyle y(t) = \frac{b x_0}{c - a} e^{-at} + \left( y_0 - \frac{b x_0}{c - a} e^{-c t} \right) \end{array} \right. \tag{4.5}\]

\[ \left\{ \begin{array}{l} x'(t)=ax(t) + by(t) \\ y'(t)=cx(t)+dy(t) \end{array} \right. \; \text{sol:} \; \left\{ \begin{array}{l} x(t) = \frac{(\lambda1 - d)x_0 + by_0}{\lambda_1 - \lambda_2}e^{\lambda_1 t} + \frac{(d - \lambda_2)x_0 - by_0}{\lambda_1 - \lambda_2} e^{\lambda_2 t} \\ y(t) = \frac{c x_0 + (\lambda_1 - a)y_0}{\lambda_1 - \lambda_2}e^{\lambda_1 t} + \frac{-c x_0 - (a - \lambda_2)y_0}{\lambda_1 - \lambda_2} e^{\lambda_2 t} \end{array} \right. \tag{4.6}\]

\[ \qquad \qquad \text{in which} \; \lambda_{1,2} = \frac{a + d \pm \sqrt{{(a+d)}^2 - 4(ad -bc)}}{2} \]

The use of Equation 4.5 requires that $c \neq a$, and of Equation 4.6 that $\lambda_1 \neq \lambda_2$. Systems of 3 or more (general case: $n$) linear ordinary differential equations with constant coefficients are solved analogously. The solutions are always sums of 3 (general case: $n$) $e$-powers. As you see, the $\lambda$’s and coefficients for $n=2$ are already quite complicated. For larger $n$, these solutions look even more complicated. However, they can still be solved.

Here ends the overview of common ordinary differential equations and their solutions. They are cookbook recipes, because we haven’t discussed how these solutions can be found. For the common user the existing solution methods are only of limited use. It costs time to learn them and they do not always work.

Checking solutions

Which does not mean that you, as a user, should believe everything a specialist tells you about them. Because, finding a solution may be difficult, but checking one is really easy. You only have to be able to differentiate. It is strongly recommended to always check solutions, also when applying the recipes above.

We will, therefore, again show the procedure for checking a solution $f$ of an ordinary differential equation, here with boundary condition $f(0)=y_0$

Checking the solution of an ordinary differential equations

Substitute $f$ on the left hand side
Substitute $f$ on the right hand side, if applicable
Check whether 1 and 2 are equal
Check the boundary condition, i.e whether $f(0) = y_0$

Step 4 is often the fastest (and easiest) way to recognize invalid solutions.

4.4 Steady states

It is often possible to deduce useful conclusions from a differential equation without having the solution itself. One of these conclusions concerns the existence and properties of steady states. The idea can be demonstrated by way of an example of a first order differential equation. The equation gives the relation between $y(t)$ and $y'(t)$. One can easily draw a graph of this relation, since the relation is given, and is not influenced by the way in which $y$ and $y'$ each depend on $t$. By analyzing the graph it is possible to draw conclusions using the normal methods of function analysis.

Example 4.6 Suppose that the length of an animal is given by the ordinary differential equation $y'(t) = 2 - 0.4y(t)$.

Plot $y'$ against $y$. The resulting figure shows how the growth rate ($y'$) depends on the length ($y$) of the animal. The process is represented by a point running in time along the horizontal axis, and the rate at which this happens can be read from the vertical axis. The arrows give an extra indication of this process, where long and short arrows indicate fast and slow change. Furthermore, when the arrow points to the right, the derivative is positive and the animal grows. When it points to the left, the derivative is negative and the animal shrinks. In this way we see that:

the maximal growth rate is found at length 0, and that it decreases with increasing $y$
growth arrests at $y=5$
when $y<5$ the animal grows and when $y>5$ the animal shrinks

Whether this model yields realistic conclusions about animal growth is a different story.

As we saw before, this ordinary differential equation is solvable, so that these conclusions could also have been drawn from the solution. The point we made here is that we can conclude from the equation itself —with much less effort— that an animal that starts at a size $<5$ will grow at a rate that gradually decreases and will stop growing at a length of 5.

Example 4.7 (Logistic growth) Take the logistic equation () that is used to describe growth of populations: $\displaystyle{y'(t)=ay(t) \left\{ 1 - \frac{y(t)}{b} \right\}}$ The relation between $y'$ en $y$ is a second order polynomial, and hence yields a parabola in the graph. Closer inspection shows that growth equals 0 at two values of $y$, namely $y=0$ and $y=b$. Furthermore, the growth rate reaches a maximum exactly half-way between these two points, at $y=\frac{1}{2} b$.

Stability

Using the last example, we will discuss an important phenomenon called stability. Consider the point $y = b$. The growth rate of the population, $y'(t)$, equals 0 there, so that the population size does not change anymore and will remain at that value once it has reached it. Suppose that, by a factor external of the growth process, a perturbation of the population occurs. A number of animals is added, or a few disappear. What will happen then? If $y$ becomes less than $b$, then (you can read it off the graph), the growth rate will be positive. Hence, the population will grow again, and thus returns to the value $b$. If the disturbance is positive ($y$ becomes greater than $b$), then $y'$ negative, and the population will drop until $b$ is reached again. In both cases, $y$ returns to the value $b$. The value $y = b$ is thus stable. A bit more precise: around that value, there is an area (which may be large or small) so that $y$ will run from any point in that area to $b$.

Hence, the technical term “stability”, is concerned with the resistance against perturbations. This is more than just a stand-still or an equilibrium. In colloquial language the difference is not made that clearly. When saying that “the patient was stable” or “the weather is stable”, we mean the stand-still of a variable process, but in “a stable traction” or “an unstable personality”, the term is perhaps more related to what we discussed here.

4.5 Simulation or numerical solution of differential equations

Simulation means imitation, and that is exactly what you do with a simulation: you imitate the process. Nowadays this is almost exclusively done using computer programs.

Euler’s method

It is often thought that the numerical solution of differential equations is only possible with specialized software. That is not the case, because it can also be done using a spreadsheet program such as Excel. Whether this solution is satisfactory depends on the type of differential equations and the accuracy you need, but using a spreadsheet program you can, using a simple method, often obtain an acceptable solution. The absolute easiest method to simulate differential equations is named after the mathematician Leonhard Euler, the same whose name is used for $e=2.71\ldots$, Euler’s number. We illustrate this method here with the differential equation which describes the decrease of artemisinine:

\[ q'(t)=-k \cdot q \tag{4.7}\]

Before we continue, we must know a boundary condition, and that in this case the value of $q$ at time $t = 0$. Suppose that value is equal to . We also need to know the value of the rate parameter $k$, and we set it at $k = 0.9 \, h^{-1}$ (the unit of $k$ is $h^{- 1}$ or $1/hour$). We will solve this differential equation numerically, but we also know the exact, `analytical’ solution, so we will be able to check how well the numerical solution approximates the exact solution. The analytical solution is

\[ q(t) = 1 \cdot e^{-0.9 \cdot t} \]

We will now apply Euler’s method to solve this differential equation numerically. To understand this method, we have to realize what Equation 4.7 means. It is a recipe to calculate the slope of the function $q(t)$ at a given time $T$, and of course that can only be done if the value of $q(t)$ is known. The problem is, of course, that the value $q(t)$ is not known in advance for an arbitrary time point $t$. However, we have a little stub to anchor our solution on. It follows from the boundary condition which says that for $t=0$, $q(t)$ is equal to 1 mM. So we know that $q(0)=1$. And we can now, using the differential equation Equation 4.7 calculate the slope at $t = 0$. It is equal to $-0.9 \cdot 1 = -0.9 \;\text{mM}/\text{hr}$. This means that at time $t = 0$, 0.9 mM artemisinine disappears per hour. We can display this knowledge graphically as shown in Figure 4.2}.

Figure 4.2: (Left) Numerical solution of the differential equation for $q(t)$ starts with the estimation of the value $q(h)$ using the boundary condition $q(0)$ and the slope in the point $(0,q(0))$. The fact that this is an estimated value is indicated by the hat on top of the $q$, $\widehat{q}(h)$. (Right) The second point $q(2h)$ is estimated by proceeding from the point $\widehat{q}(h)$. The fat line indicates the true, but unknown function $q(t)$, and the coloured lines are the tangents with the calculated slopes.

It is now probably already clear how we can proceed from there. Since, from this figure you see that, using the initial value $q(0)$ and the slope $q'(0)$, we can estimate the value $q(h)$, where $h$ is a small step in time, like $h=0.01$. An estimate of $q(h)$ is

\[ q(h) \approx \widehat{q}(h) = \underbrace{q(0)}_{\textrm{start}} \quad + \quad \underbrace{q'(0)}_{\textrm{slope}} \quad \times \underbrace{h}_{\textrm{step size}} \]

If we substitute the values for $q(0)$, $q'(0)$ and $h$, we obtain $\widehat{q}(0.01) = 0.991$. The fact that this is an estimate is indicated by the hat $\widehat{q}$. For an arbitrary time $T$ we are able to estimate the value of $ q(T+h)$ if we know the value of $q(t)$:

\[ q(T+h) \approx \widehat{q}(T+h) = \widehat{q}(T) + \widehat{q'}(T) \cdot h \]

Hence, now that we have an estimate of $q(0.01)$, we can also make an estimate a value $q(0.01 + h)$ a little bit later. If we take the same time-step $h = 0.01$ again, this will be $q(0.02) \approx \widehat{q}(0.02) = \widehat{q}(0.01) + \widehat{q'}(0.01) \cdot 0.01$. The derivative $q'(0.01)$ can be estimated as $q'(0.01) \approx \widehat{q'}(0.01) = \widehat{q}(0.01) \cdot -0.9 = 0.991 \cdot -0.9 = -0.8919$. The estimated value for $q(0.02)$ will then be $\widehat{q}(0.02) = 0.982081$.

In this way we can calculate a chain of estimates for $q(t)$ for the time points 0, 0.01, 0.02, etc. The estimate at a specific time point is always calculated from the estimated values for $q$ and $q'$ of the previous time point. This estimate of the solution of the ordinary differential equation is not a function definition, but a table with three calculated columns. The first column contains the time points, the second the estimated values of $q$, and the third the estimated value of $q'$. Such a table is shown below, including the formulas. This result is called a numerical solution.

Formula for time	Time	Formula for estimated value	Estimated value of $q$	Formula for estimated slope	Estimated slope, $q'$
$0$	0	$q(0)$	1	$q'(0) = k \cdot q(0)$	-0.9
$h$	0.01	$\widehat{q}(h) = q(0) + h \cdot q'(0)$	0.99100	$\widehat{q'}(h) = k \cdot \widehat{q}(h)$	-0.89190
$2h$	0.02	$\widehat{q}(2h) = \widehat{q}(h) + h \cdot \widehat{q'}(h)$	0.98208	$\widehat{q'}(2h) = k \cdot \widehat{q}(2h)$	-0.88387
$3h$	0.03	$\widehat{q}(3h) = \widehat{q}(2h) + h \cdot \widehat{q'}(2h)$	0.97324	$\widehat{q'}(3h) = k \cdot \widehat{q}(3h)$	$\ldots$
$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
$800 h$	8	$\ldots$	0.00072	$\ldots$	$\ldots$

Mathematically, this procedure is not quite right: the numerical solution is also explicitly an approach of the mathematical solution. Because, during a time step $h$ we pretend that the growth rate is constant, and is equal to the rate at which it was at the beginning of the corresponding time step. That’s only approximately correct, because also during that short period $h$, $q'(t)$ will change! The error can be kept small by keeping the time step $h$ small. However, for many differential equation, the larger the simulation, the greater the deviation. Therefore, it is always a good idea to test a few values for $h$ to investigate the influence of this parameter on the numerical solution. Ideally, this influence is negligibly small. Euler’s method is a so-called first order approximation, because we use as a straight line in the estimate to approximate the change of the function $q$ over the time-fragment $h$. There are other method that use second- and higher-order approaches to approximate the function. Furthermore, the size of the time steps may also be varied, for instance by making it small in time regions where the solution changes quickly, and larger when it changes only little. Using these methods, the numerical approach can be greatly improved and also be made faster. These methods are of course more complex than the Euler’s method. Software for solving differential equation usually contains multiple algorithms to approximate a solution for different types of differential equation. Examples of such software are Mathematica (www.wolfram.com/mathematica, commercial, academic licenses available), Maxima (maxima.sourceforge.net, free, open source) and Scilab (www.scilab.org, free, open source). By the way, the range of applications of this software is much broader than just for solving differential equation.

From the description above, an essential property of computer simulation is clear: to calculate the time course of the solution, the solution itself is not needed. The ordinary differential equation itself suffices. That means you can obtain a solution for any ordinary differential equation. Given that many ordinary differential equations can not be solved analytically, this is a valuable conclusion.

There are also disadvantages: one only obtains a graph of the solution, not an analytical equation. Also, simulation can only give a particular solution, and not the general solution. Both points lead to a loss of transparency, and to learn how the solution changes with different parameter values, for example, multiple simulations have to be run.

4.6 The SIR model for spread of infectious disease

The mathematical modeling of spreading of infectious disease like common cold, COVID-19 or Ebola is an important part of epidemiology. Such models are, for example, used to monitor and predict the effects of partial immunity, vaccination programs and other interventions. Here we will use a simple and very effective model for the spreading of infectious disease called the SIR model. This model, and variants of it are used to monitor epidemics like COVID-19. It consists of a set of three differential equations for three variables called $s$, $i$ and $r$. These three variables indicate three different fractions of persons in a population that are in either of three conditions or groups, namely

Susceptible: Individuals in the Susceptible condition have never been infected before, and they are able to catch the disease. Once they have it, they go into the Infected condition.
Infected: Individuals in the Infected condition can spread the disease to individuals in the Susceptible condition. The time they spend in the Infected state is the infectious period, after which they enter the Recovered condition.
Recovered: Individuals in the Recovered condition are assumed to be immune for life (or for the duration of the epidemic). So, they stay in this condition.

Every person in a population is at any time in one and only one of these conditions but can transit from one group or condition to another. Let’s start by writing down these transitions as if they were chemical reactions, concerning S-type, I-type and R-type molecules (Figure 4.3).

Figure 4.3: Transitions in the SIR model. Dots indicate persons. The population is split virtually and a person belongs to one of the classes or conditions *Susceptible* (S), *Infected* (I), or *Recovered* (R).

This transition scheme is not only a graphical display of the transitions, but also a reminder of the fact that we can choose to model transitions as if we are working with mass-action kinetics, in this case of truly irreversible transitions. Ideally, the population is an ensemble of independently reacting “particles”. For example, we can model the rate of change of persons going from Infected to Recovered condition as a “reaction” of which the rate only depends on the fraction $i$ of I-individuals:

\[ v_{\text{recovery}} = \gamma \cdot i \]

The rate of infection, i.e. rate of persons transiting from Susceptible to Infected condition is dependent on the number of Susceptible persons, but also on the number of Infected persons. We can model the rate of the transition $S \longrightarrow I$ from Susceptible to Infected as:

\[ v_{\text{infection}} = \beta \cdot s \cdot i \]

We now have differential equations for $s$ and $r$, and we still need one for $i$. The rate of change of $i$ is, of course:

\[ \dd{s}{t} = -v_{\text{infection}} = - \beta \cdot s \cdot i \tag{4.8}\]

That with which $I$ changes is \[ \dd{i}{t} = v_{\text{infection}} - v_{\text{recovery}} = \beta \cdot s \cdot i - \gamma \cdot i \tag{4.9}\]

and finally, that with which $r$ changes: \[ \dd{r}{t} = v_{\text{recovery}} = \gamma \cdot i \tag{4.10}\]

The SIR model consists of these three equations. These are called coupled differential equations because the variables occur in each other’s differential equations. In a coupled system of differential equations, finding a solution for one of the variables often means that you have also found a solution for the others. The model is characterized by only two parameters, namely $\beta$, a constant indicating the effectiveness of the infection process when an infected individual encounters a susceptible individual, and the constant $\gamma$, which determines the rate of recovery of infected individuals. Additionally, the boundary- or initial conditions are important. The number of resistant individuals in the population at the start of the epidemic can be changed by varying the initial value of the variable $r$, for example. In this way, you can study the effect of immunity in the population on the time-course of the epidemic.

The reproduction numbers $R_\text{eff}$ and $R_0$

Taking a closer look at Equation 4.9 we see that it can also be written as

\[ \dd{i}{t} = i {\left( \beta \cdot s - \gamma \right)} \]

which immediately shows that $i$ will increase if $\beta \cdot s - \gamma > 0$ and decrease otherwise. This condition can be re-written as

\[ \text{if }\frac{\beta}{\gamma} s > 1\text{ then }i\text{ decreases} \]

The number $\frac{\beta}{\gamma} s$ is called the effective reproduction number $R_\text{eff}$ and $\frac{\beta}{\gamma}$, i.e. when the entire population is susceptible and $s=1$, is called the basal reproduction number $R_0$. In principle, in the SIR model, $\beta$ and $\gamma$ are constants, but in practice the efficiency of infection $\beta$ changes (decreases) when measures are taken to decrease the spreading of an infectious disease. Therefore, the effective reproduction number $R_\text{eff}$ is an important quantity for monitoring an epidemic. See for example the estimated $R_\text{eff}$ of COVID-19 in the Netherlands in Figure 4.4.

Figure 4.4: The estimated effective reproduction number during the first half year of the COVID-19 epidemic in the Netherlands. Social restrictions were imposed in the first half of March. Source: RIVM.

Simulation of an epidemic

An algebraic solution for the system of differential equations of the SIR model is difficult to find, but we can easily solve it numerically by the method of Euler, described above. You will have to determine a time-step $\Delta t$, and have to make a table with 7 columns: time, $s$, $i$, $r$, $\Delta t \cdot \tdd{s}{t}$, $\Delta t \cdot \tdd{i}{t}$, and $\Delta t \cdot \tdd{r}{t}$, as in Figure 4.5.

Figure 4.5: Spreadsheet file for the simulation of an infection using the SIR model. The free parameters are shown in yellow. The formula in the gray-green cell is shown in the *formula-bar*.

4.7 Exercises

Exercise 1.

The growth of an organism is given by the Von Bertalanffy growth model. The volume $V$ of the organism is a function of the age $a$. $\alpha > 0$ and $\beta > 0$ are two parameters. The growth is described by the following differential equation: \[ V'(a) = \alpha {V(a)}^{\frac{2}{3}} - \beta V(a) \]

The first part shows that the intake of food is proportional to the surface (surface is a square of the distance and volume a cube, so surface is volume to the power $\frac{2}{3}$) and the second part shows that the maintenance is proportional to the volume.

a. Calculate the equilibrium ($V^*$)

Tip

Change variables to obtain growth $V'$ as a function of size $V$, e.g. $g = V’(a)$ and $s = V(a)$ and then draw the graph for $g(s)$.

b. Calculate the volume where maximal growth is reached ($V'_{max}$)

c. Sketch the graph of $V'$ against $V$ for $\alpha =1$ and $\beta =1$. Is the equilibrium stable?

Exercise 2.

A patient gets medication via an infusion. The amount of medicine in the patients body is called $m$. The patient gets $a$ mg medication per hour. The velocity of the clearance of the medicine can be the result of different processes, but is usually a first order process (here indicated by $b$ in mg per hour). The amount of medicine in the patients body ($m(t)$) obeys the following ODE: \[ m'(t)= a - b \cdot m(t) \]

The amount of medicine at $t=0$ is $m_0$.

a. Give the solution of this ODE (the function $m(t)$).

b. Assume that the clearance is 0.1 mg per hour and the amount of medicine at the beginning is 2.5 $\text{mg}$. How much medication should you add per hour to be sure that the amount of medicine in the patients body is always at least 2 $\text{mg}$?

Exercise 3.

A substrate (S) is converted into an intermediate (X), and subsequently into a product (P) with mass-action kinetics. The substrate is fixed at 100 $\text{mM}$ and the product is directly converted (fixed at 0 $\text{mM}$).

\[ S \rightleftharpoons X \rightleftharpoons P \]

Hint: Mass action kinetics means that the rate is proportional to the concentration. For example if S is converted to X the forward rate is $k_1 \cdot [S]$ and the backward rate is $k_{-1} \cdot [X]$. In this question we have an intermediate, thus there are 2 such reactions! The $k$’s are constants and you can leave them in you answer (since you don’t know the value).

a. Give the differential equation for the change in X, $\dd{(X)}{t}$. For this you have to write down with which net rates $X$ is produced and consumed. Finish by filling in $S$ and $P$.

b. Solve the differential equation to find the concentration of X in time with $X_0=0$.

c. How long does it take for X to reach half of it’s steady state value?

Exercise 4.

A researcher studies a transcription factor in a fast growing human cell line. She concludes that the kinetics of the factor can be described by the following differential equation: \[ \frac{\mathrm{d}x}{\mathrm{d}t} = a \left( 1 - \frac{x^n}{b^n + x^n} \right) - k x \] where $x$ is the amount of transcription factor in the cell, $t$ is time, and $a$, $b$, $n$ and $k$ are positive constants.

a. Indicate which parts of the formula symbolize the synthesis rate of the transcription factor and which parts the rate of decrease.

b. What are the maximal and minimal rates of synthesis?

c. Sketch the rate of synthesis as a function of the amount of transcription factor when $n \gg 1$

d. By which molecular mechanism could the value of $n$ be determined?

e. What is the name of this kind of regulation?

f. In the graph made above, indicate the steady states. Also indicate which of these steady states are stable and which are unstable.

Exercise 5.

Logistic growth is often used for modeling population growth. It is described by the following differential equation:

\[ \dd{X}{t} = X \left( 1 - \frac{X}{100} \right), \qquad X \geq 0 \]

a. Sketch the graph of $X'(t)$ against $X(t)$. To do that, find the intersections with the axes and the extrema.

b. If we start with a population size of 1, what will be the eventual population size? And what if we start with a population of size 1000?

c. Which of the following growth equation describe logistic growth? , which of equations are a solution to the logistic growth equation?

$X(t) = 100 \cdot e^{t/2}$
$X(t) = \frac{1}{1+e^{-t}}$
$X(t) = \frac{100}{1+e^{-t}}$

Exercise 6.

We have the following differential equation: \[ x'(t) = \frac{5 x(t)^2}{4 + x(t)^2} - x(t) \]

a. Sketch the graph of $x'(t)$ against $x(t)$. To do that, find the intersections with the $x(t)$ axis.

b. What will be the eventual value of $x$ when $x(t)$ starts at $x=\frac{1}{2}$, $x=3$, or $x=1$? Explain why the starting point $x=1$ is special.

c. Somebody asserts that the following function is a solution to the differential equation. Is this assertion correct? \[ x(t) = \frac{1}{1 + e^{-t}} \]

Metabolic networks

Exercise 7.

Consider the following network (reactions are defined from left to right, but are reversible for the kinetics):

Figure 4.6: Reaction network. The kinetics (mass-action kinetics) are reversible, but the reactions are defined from left to right.

a. Write down the differential equations of $B$, $C$ and $D$ in terms of the rates ($v_1$, …, $v_5$).

b. Assume there is steady state for the intracellular metabolites $B$, $C$ and $D$ (all differential equations are $0$). Which rates are equal?

c. Assume mass action kinetics with $A = 100$ and $E = 0$. Write down the differential equation in terms of the metabolite concentrations ($B$, $C$, and $D$) and the kinetic constants ($k_1, k_{-1}, \dots$).

Exercise 8.

Consider the following network:

Figure 4.7: Reaction network. The kinetics (Michaelis-Menten kinetics) are irreversible, and the reactions are defined from left to right.

a. Write down the differential equation for $B$ in terms of the rates ($v_1$ and $v_2$).

b. Assume that the rates depend hyperbolically on their substrates and not at all on their products: \[ v_1 = \frac{v_{\max 1} [A]}{[A] + K_1}, \quad v_2 = \frac{v_{\max 2} [B]}{[B] + K_2}. \] Moreover, assume that $[A] = 2 \,\text{M}$ and $[C] = 0 \,\text{M}$ and that $v_{\max 2} > v_{\max 1}$. Sketch $v_1$ and $v_2$ as a function of $[B]$.

c. Sketch $v_1 - v_2$ as a function of $[B]$. How can you see the steady state of $[B]$ in this graph and for what parameter values does this steady state exist?

d. Sketch $B$ versus time for

\[ v_{\max 1} = 4, \quad K_1 = 1, \quad v_{\max 2} = 5, \quad K_2 = 1, \quad [A] = 2, \quad B(0) = 0. \]

The SIR model

Exercise 9.

This question concerns the COVID-19 epidemic.

a. Estimate the basal reproduction number $R_0 =\frac{\beta}{\gamma}$ from Figure 4.4 under the condition where no social restrictions are imposed.

b. The COVID-19 virus will not spread in the population, or better: it will die out, if $R_\text{eff} < 1$. What is the fraction of the population that needs to be resistant to obtain an $R_\text{eff}$ that stops a COVID-19 infection from spreading under the condition that no social restrictions are imposed, in society as it was before March 2020 with the same value for $\beta$? ().

Exercise 10.

In the popular press, the reproduction numbers $R_0$ and $R_\text{eff}$ are often said to be the average number of new infections that an infected person causes. Here, we want to check this assertion. Clearly, the number of infections caused by an infected person depends on the period of time that this person is ill. Let’s call the probability that a person is still ill at a certain time $t$ after infection $P_{i}(t)$, to distinguish this from $i(t)$, the fraction of infected persons in a population. The function $P_{i}(t)$ is subject to the following differential equation: \[ \dd{P_{i}}{t} = -v_\text{recovery, person} = -\gamma P_{i} \tag{4.11}\]

i.e. the probability that a person recovers is proportional to the probability that the person is still ill. (You should see the fact that this equation concerning a single person implies the rate equation for $v_\text{recovery}$ valid for the fraction of the population $i$ that recovers per time unit).

a. Solve Equation 4.11 to describe $P_{i}(t)$ explicitly as a function of time.

b. Using the previous solution, calculate the avarage time that a person is ill.

c. The rate at which an infected person infects others equals $s \beta P_{i}(t)$. Calculate the average number of persons infected by an infected person (assuming $s$ hardly changes during the time that a single person recovers).

Exercise 11.

Set up a simulation in a spreadsheet program of the SIR model using the parameters as given in Figure 4.5. Again, we will study the COVID-19 data.

a. Assume that the recovery of a COVID-19 infected person takes on average two weeks. During this period the person is infectious. Estimate $\gamma$ using “weeks” as a time unit, and calculate $\beta$ from this number and the number $R_0$. Substitute both numbers in the spreadsheet.

b. Estimate the peak fraction of infected persons, and given that the total population of the Netherlands is 17 million, the peak fraction of infected persons. If 1% of infected persons need hospitalization, what would have been the peak number in those persons (note that we do not take into account the fact that on average hospitalized persons needed much longer than two weeks to recover).

We have a solution for a mathematical equation if substitution of that solution yields the same expression on both sides of the equation. Is $x=2$ a solution to $x^2=5x$? The left hand side yields 4 and the right hand side 10, so no!. $x=5$ is, however, a solution.↩︎

Formula for time	Time	Formula for estimated value	Estimated value of \(q\)	Formula for estimated slope	Estimated slope, \(q'\)
\(0\)	0	\(q(0)\)	1	\(q'(0) = k \cdot q(0)\)	-0.9
\(h\)	0.01	\(\widehat{q}(h) = q(0) + h \cdot q'(0)\)	0.99100	\(\widehat{q'}(h) = k \cdot \widehat{q}(h)\)	-0.89190
\(2h\)	0.02	\(\widehat{q}(2h) = \widehat{q}(h) + h \cdot \widehat{q'}(h)\)	0.98208	\(\widehat{q'}(2h) = k \cdot \widehat{q}(2h)\)	-0.88387
\(3h\)	0.03	\(\widehat{q}(3h) = \widehat{q}(2h) + h \cdot \widehat{q'}(2h)\)	0.97324	\(\widehat{q'}(3h) = k \cdot \widehat{q}(3h)\)	\(\ldots\)
\(\ldots\)	\(\ldots\)	\(\ldots\)	\(\ldots\)	\(\ldots\)	\(\ldots\)
\(800 h\)	8	\(\ldots\)	0.00072	\(\ldots\)	\(\ldots\)