# Differential Equations

## An interactive tutorial

*In this chapter, you’ll learn how to model systems using ***ordinary differential equations*** and ***Python. ***Along the way, you’ll learn two famous systems of differential equations: the ***predator-prey model*** of ecosystems, and the ***Forbes model of metabolism*** (you might be surprised how your own body works!).*

## Warm-up: the predator-prey model

Ahh ... another lunchtime out in nature. Just you, your binoculars, and a picnic lunch. Bird-watching truly is bliss. But wait ... I hear a bird rustling under that bush! What species is it?

Be careful! Advers are known to bite!

Aww, so it is! The first one today.

Over the past couple of years, you’ve kept diligent records on these two local bird species. Here’s the result of your tallies in chart form:

You think? I’m a bit worried ... aren’t the Midgens dying out? Will they be completely replaced by the Advers?

Indeed! Are the Midgens dying out? Will they be completely replaced by the Advers?

Your records are not long enough to see the long-term trends ... what you need is a kind of **mathematical model** to fit the data to.

However did you guess? Yes, differential equations can model all sorts of things, and are especially great when you don’t have much data (as is the case here).

So let’s start modelling the avian ecosystem from first principles: the birds and the bees. Left to themselves, Midgens beget Midgens. Those baby Midglets are adorable. 🐣 We can express this wonderful fact of life in our first **Ordinary Differential Equation, **or **ODE:**

Isn’t it fun? This is not a course in calculus, but here’s the one thing you need to know: the letter $d$ can be read as “change in”. For example, we’re using $🐦$ to stand for the number of Midgens, so $d🐦$ can be read as “change in the number of Midgens”.

Given that $t$ stands for time, what does $dtd🐦 $ represent?

No. It’s true that “speed” is one way to talk about “change rate”: the rate of change in position. But we’re instead talking about the rate of change of the Midgen population, not the position of a Midgen. $d🐦$ is change in Midgen count, and $dt$ is change in time. Divide one by the other, and you get change in Midgen count over time.

So the equation above, in English, can be read as “The change in the number of Midgens over time is proportional to how many Midgens there are.” Here, $α$ is the proportionality constant, and represents the natural growth rate of the Midgens.

This is perhaps the simplest ODE, but it’s already interesting. For example, what would happen if $α=0$?

If $α=0$, our equation simplifies to $dtd🐦 =0$. In other words, the change over time is zero, so the population stays at its initial population.

Indeed, the equation simplifies to $dtd🐦 =0$, meaning there is no change over time.

Let’s say we set $α=0.002$. That is, each day, each Midgen begets $0.002$ new Midglets. Or to put it in more familiar terms, each Midgen has a baby every $500$ days.

Run this model in your head. How does the population grow over time?

Luckily, we don’t have to run it in our head! Instead, we can run it with an “ODE solver” like `scipy.integrate`

. Later, you’ll learn how to do this; for now, here’s the graph:

Well, look at that. Exponential growth! Without any competition, the Midgens will multiply unchecked.

## A system of two equations

Indeed. This is not an accurate model. So far you’ve learned how to model change in one variable (🐦), but the world is more complex. Let’s learn how to deal with a second variable: the number of Advers (🦅).

The Midgens are a prey species, kept in check by the predator Advers. They eat those innocent Midglets. 🍳 We must add this second fact of life to our mathematical model. Here’s the full model, called the **predator-prey model:**

You can interpret this with some thought.

The term $🐦🦅$ is the number of Midgens multiplied by the number of Advers. You can think of this as the number of interactions between the two species. $β$ says how frequently those interactions result in a dead Midgen, and $δ$ says how much those dead Midgens are converted into new Advers.

Now, think about the term $−γ🦅$. It’s used in the second equation, describing the change in Advers over time ($dtd🦅 =…−γ🦅$). What could that term $−γ🦅$ represent?

Right, it could include birth rates and death rates. The full term is $−γ🦅$, representing processes that change the number of Advers in proportion to the population size. Births and deaths are both proportional to the population.

I don’t think so: if this represented the number of Advers killed by Midgens, it would include a $🐦$. But the full term is $−γ🦅$, representing a process that changes the number of Advers in proportion to just the current population size. Births and deaths (due to old age) are both proportional to the population size, so this term includes births and deaths.

To actually use this model, you need to estimate these parameters. But luckily, you have also studied the birth, death, and kill rates, and you use these to come up with:

(Note that $γ$ is positive! A population of Advers on their own will die off, rather than grow.)

Now, try again to run this model in your head. It’s pretty difficult, isn’t it?

Okay, here’s your graph:

The population change is cyclical! You’ve now seen two kinds of system: one that never reaches a stable state, and another that reaches a kind of cyclic stability.

Finally, your observations make sense. You saw a dip in the Midgen group and a growth in the Advers, but soon the trends will reverse again: the Advers will become victims of their own success, destroying their food supply and killing each other. The circle of life continues.

## Running ODEs in Python 🐍

Oh yes! We can translate a model into Python by defining an initial state, and a function that describes the rate of change. Here’s a Python equivalent of our first model, which grows exponentially:

```
s0 = 100 # Starting number of Midgens
def rate_of_change(s, t):
return 0.002 * s
```

We can then use the `scipy.integrate`

package to run this model:

```
import numpy as np
from scipy.integrate import odeint
t = np.linspace(0, 365*3)
timeline = odeint(rate_of_change, s0, t)
```

Finally, we can graph the timeline with:

`plt.plot(days, timeline)`

That’s basically all the Python you need to know! For the complete runnable code, check out our Colab notebook.

## The Forbes model of metabolism

*And now for something completely different.* You’ve learned the basics of ODEs and how to run them. Now you’ll use them to model a very real-world complex system: **your own body! 🚶♀️**

This was a preview of chapter 2 of Everyday Data Science. To read all chapters, buy the course for $29. Yours forever.

Last year, my bathroom scales said I weighed $55$ kg. A ruler says I am 1.82 meters tall. This put my BMI at $17$, which is classified as underweight. I set myself the target of reaching at least $63$ kg, bringing my BMI into the “healthy” range. But how could I achieve this?

We’ve all been told: If you eat more calories than you burn, you will ...

No, it’s the other way around.

And this advice is basically true, but makes it sound simple. We often talk about the $2000$ kcal diet as if it were a universal truth, but each person is different! How many calories are *you* burning right now?

Again, you need a model. Given a target weight, you can calculate your *personal* caloric intake using a **‘weight change’ ODE.**

We’ll call the amount you intake $EI$, and the amount you expend $EE$. The difference $EB=EI−EE$ is called your **energy balance.** If your intake equals your expenditure, $EB=0$, you’re perfectly balanced, and your weight will be stable.

We’ll measure all of these in **kilocalories **(usually just pronounced “calories”) per day. Yesterday, I ate $2500$ kcal, and I burned $2700$ kcal. What was my energy balance for the day?

Not quite, it’s the other way around: $EB=EI−EE=2500−2700=−200$ kcal.

Not quite: it’s $−200$ kilo*calories*, not kilo*grams*.

So these variables are measurements of *energy*, but ultimately we’re concerned with *weight*. How are energy and weight related? Einstein had important opinions here ⚛️, but instead we’ll use a model developed by Gilbert Forbes. Forbes says, if you squint, your body (measured in **kg**) is made up of two kinds of tissue :

Now, check the nutritional fact label 🏷️ on the back of your neck:

```
HUMAN BODY (approx. 50 servings)
Fat tissue ............. 9400 kcal/kg
Lean tissue ............ 1800 kcal/kg
```

It’s been there since birth! (Your figures might differ a little bit from the Forbes model, but probably not by much.)

This label is helpful alongside a modern bathroom scale, which reports weight and percentage fat. Let’s try that out. Kevin’s scale says he weighs $110$ kg, of which $10$ kg is fat. What is Kevin’s total calorific content?

The correct calculation is:

(Incidentally, this means you could use Kevin to heat one liter of water to around $274000°C$. But just because you *could*, doesn’t mean you *should*.)

Now, when you eat, should your body build more Fat tissue or more Lean tissue? We call your body’s decision $p$: it’s a *proportion* between $0$ and $1$. If $p=0$, your body only builds/burns Fatty tissue. If $p=1$, your body only changes Lean tissue.

But how does your body make this decision, $p$? Forbes says it’s based on how much fat you have. The model looks like this:

$C=54.3$ is a magic constant. Let’s sanity-check this function. Say you have no fat on you, and then you eat a meal. Does your body use the meal to build Fat tissue or Lean tissue?

Actually, if you have no fat, $p=C+FatC =C+0C =1.$ Check above: $p=1$ means your body only builds lean tissue.

On the other hand, as $Fat$ tends to $∞$, your body only burns/builds Fat tissue.

We can now write two ODEs describing how your body changes:

Here, your net energy $EB$ is being shared between your Lean and Fat tissue, and then converted from kcal to kg.

Now recall the original equation, $EB=EI−EE$. The value $EI$ is just how much we decide to eat, and we will play with $EI$ to obtain our target output weight.

But calculating the energy expenditure $EE$ is more complicated. Let’s simplify the problem by assuming that you never move.

No, you can keep your shape. Bed-bound, your body burns fuel just to stay alive. This is called your **basal metabolic rate****,** or $BMR$, measured in **kcal per day**. Here is Forbes’ model:

Don’t panic! 🙂 The first line here just says it takes energy to *maintain* current tissue. The second line says it takes energy to *build* new tissue.

Now just consider the last term, $BMR=…+0.14EI$. How would you interpret this term?

That’s not it. (That might be written as something like $EB=0.14EI−EE$.) The correct interpretation is that increasing your intake $EI$ (by consuming food) increases your metabolic rate $BMR$ (i.e., burns energy).

Yes, it even takes energy to consume energy! (This is called the **Thermic Effect of Food**.)

Now you should understand why bears store fat when hibernating: to store the same energy as lean tissue, they would need *five* times as much space! And to maintain that lean tissue would cost almost *seven* times as much! 🐻❄️

Finally, to get your energy expenditure $EE$, Forbes says you must multiply your basal metabolic rate by your **physical activity level, **called $A$. Now, how many hours do you exercise per day?

Impressive. Your $A$ is between $2$ and $2.4$.

Okay, your $A$ is between $1.7$ and $2$.

Okay, your $A$ is between $1.4$ and $1.7$. Try to walk around a bit.

The activity level is an abstract number between $1.4$ (no exercise) and $2.4$ (multiple hours of exercise every day).

Note that “no exercise” is still above $1$: you probably move around more than what’s considered truly “basal”.

Note also that exercise is not just *added* to your metabolic work — it *multiplies* your metabolic work! By exercising, you make your body spend more even when it’s resting!

Finally, we solve for the energy balance $EB$. I’ll spare you the algebra! Here’s the ugly but functional equation:

We now have everything we need to run the model. But again, try to run the model in your head first! For example, given a high-calorie diet like $EI=4000$ kcal, my weight will start increasing, but what will happen to my weight eventually?

Well, let’s find out! Using the trusty ODE solver, here are various plots with different values of $EI$:

It turns out that weight always approaches a stable state! My personal goal is to reach $63$ kg. How much should I consume per day?

Check out the graph: at $4000$ kcal/day (the blue line), my weight would grow to over $80$ kg. It’s sufficient to eat $3000$ kcal/day (the orange line), where my weight would reach just over $63$ kg.

Check out the graph: at $2000$ kcal/day (the red line), my weight would actually drop further, to a dangerously low BMI. At $3000$ kcal/day (the orange line), my weight would reach just over $63$ kg.

Right! With $3000$ kcal/day (the orange line), my weight will reach just over $63$ kg.

I was able to follow this regimen (although it was quite challenging) and move much closer to a healthy weight. It takes a lot of effort to maintain your health, but by using all of the tools at your disposal, you can give yourself the best chance to meet your goals.

You’ve now seen three types of system: one system with no stable state, one system with cyclic behavior, and one system that settles to a stable state. Differential equations are very expressive!

**End notes:**

⭐ For the complete runnable code, check out our Colab notebook.

⭐ The predator-prey model was originally designed by Lotka and Volterra.

⭐ They’re called *ordinary* differential equations because we only talk about change with respect to *time*. Or, more precisely, a single independent variable, which in all our models was the time $t$.

⭐ Roughly speaking, `odeint`

just repeatedly runs our `rate_of_change`

function, and sums up the outputs. Although it does clever stuff to be more accurate and efficient. You could look into Runge-Kutta methods.

⭐ For more details on the weight model, peruse “Modeling weight-loss maintenance to help prevent body weight regain” or “Quantification of the effect of energy imbalance on bodyweight”.