Modeling Life: Foundations

Herman Wouk was researching for a novel he planned to write about World War II. He interviewed several physicists at Caltech who had worked on the bomb, including Richard Feynman. After the interview, Feynman asked Wouk if he knew Calculus. “No,” said Wouk. To which Feynman replied, “You’d better learn it. It’s the language God speaks.” I came across this quote from the wonderful book Infinite Powers.

I recently read Alan Garfinkel’s wonderful book Modeling Life and watched his accompanying YouTube lectures, which I highly recommend as the best introductory book for learning calculus. Together, they teach how to apply the core concepts from Calculus and Linear Algebra to create models that can predict the behavior of dynamic systems.

In this post, I’ll scratch the surface of what goes into creating a model for a simple Shark-Tuna ecosystem. 

What would happen if you put a group of sharks and a group of tuna in a large patch of ocean? Naturally, the sharks would begin to prey on the tuna. But after some time, would the ecosystem end up with more sharks and less tuna, or less sharks and more tuna? How do we go about predicting the future state of this predator-prey system?

Condition of a system at a given time is represented by its state variables. State variable is a number that characterizes a system at a particular time. The state variables here are the number of sharks (S) and the number of tuna (T). State variables are a function of time. Let there be 100 sharks and 30 tuna at the start.

The state of the Shark-Tuna ecosystem can be represented as a 2-dimensional vector (S, T), where S represents the number of sharks and T represents the number of tuna. At the start, this vector has a value of (100, 30). The question is, where will this state vector end up after some time has passed?

To answer this question, we need to understand what factors cause the shark population (the S variable) and the tuna population (the T variable) to change over time. Let’s start by examining what influences the tuna population.

The tuna population increases when tuna reproduce and give birth. The tuna population decreases when tuna are eaten by sharks. To simplify, let’s ignore the natural death of tuna for now. The positive term affecting the tuna population is the per-capita birth rate multiplied by the number of tuna. The negative term is the rate at which sharks encounter and successfully prey on tuna.

Change in Tuna = (tuna_birth_rate * T) – (probability_shark_win * S * T)

S*T above assumes that each shark has an equal chance of encountering each tuna, based on how many of each are in the ocean.

The shark population increases when sharks successfully catch and eat tuna. The negative term in the tuna population equation now becomes the positive term for the shark population. Since sharks are much larger than tuna, let’s assume that each shark needs to catch and eat 2 tuna in order to give birth to one new shark. The shark population decreases through natural death.

Change in Shark = ((1 / tuna_per_shark) * probability_shark_win*S*T) – (shark_death_rate * S)

The two equations presented above are called a model or a change equation. They are also referred to as differential equations. The model instructs how the state variables (the shark and tuna populations) change over time.

In this model, S and T represent the state variables, as their values change over time. In contrast, the quantities tuna_birth_rate, probability_shark_win, tuna_per_shark, and shark_death_rate are constants, and are referred to as the parameters of the model. We will set the following values for the model parameters:

1. tuna_birth_rate = 0.5
2. probability_shark_win = 0.01
3. 1 / tuna_per_shark = 0.5 where tuna_per_shark value equals 2
4. shark_death_rate = 0.2

Let’s run the change equations using the parameter values we just defined and see what results we get in one run. What observations can you make from the output?

Time	S	T	S-change	T-change	S-change * delta_t	T-change * delta_t
0	100	30	-5	-15	-0.05	-0.15
1	99.95	29.85	…	…	…	…

The change in tuna population according to the tuna change equation is -15. To get the actual change over a small time interval, I multiplied this value by a delta_t of 0.01. I will explain the reasoning for using a small time interval shortly.

The resulting value of -0.15 tells us that the tuna population is decreasing. If we add this change of -0.15 to the initial tuna population of 30, we arrive at 29.85 tuna for the next time cycle. In real life, the tuna population would be in the billions, so using decimal numbers here is just for the sake of this illustrative example.

The change in the shark population according to the shark change equation is -5. Multiplying this by a delta_t of 0.01 the resulting change of shark is -0.05. If we add this change of -0.05 to the initial shark population of 100, we arrive at 99.95 sharks for the next time cycle.

The updated values for the state vector (S, T) after applying the changes are (99.95, 29.85). This updated state vector is shown in the second row of the table above. The change vector, with values (-0.15, -0.05), represents the changes in the shark and tuna populations, respectively. Change vector is pointing left and down as both the change values are negative.

Change vector exists in the change vector space. It instructs the state variables living in the state vector space to update themselves for the next cycle. The table below shows the changes in the shark and tuna populations over 4 time cycles.

Time	S	T	S-change	T-change	S-change * delta_t	T-change * delta_t
0	100	30	-5	-15	-0.05	-0.15
1	99.95	29.85	-5.07	-14.91	-0.051	-0.149
2	99.90	29.70	-5.14	-14.82	-0.051	-0.148
3	99.85	29.55	-5.22	-14.73	-0.052	-0.147
4	99.80	29.41	-5.29	-14.64	-0.053	-0.146

Looking at the trend, it appears that both the shark and tuna populations are declining. One might be tempted to extrapolate and conclude that after 10,000 cycles, their populations will continue to decrease. However, that would be an incorrect assumption.

Systems with feedback loops, like this predator-prey model, can behave in counterintuitive ways that make simple extrapolation misleading. The long-term dynamics of such systems are often more complex and cannot be easily predicted by just observing the short-term trends.

I ran this simulation for 10,000 cycles, and I was amazed to see the shark and tuna populations continue to oscillate over time. At the start, there were more sharks feeding on fewer tuna. Both populations continued to decline as a result.

However, as the shark population decreased, the pressure on the tuna population was reduced, allowing the tuna to grow in number. As the tuna population increased, the sharks found more prey available, causing the shark population to start growing again. More sharks puts more pressure on tuna population causing tuna population to decline. The cycle then repeats.

The key to understanding this system is time delays. Sharks eat tuna, so the shark population grows and the tuna population diminishes until we get to a state where there are many sharks and very few tuna. The shark growth was caused by the previously high tuna levels, but now the tuna levels are low. The delayed shark growth has created a high- shark/low-tuna state, which means that the shark population will then decline, because due to the low tuna levels, very few sharks will be born and/or survive to maturity. This shark crash then takes the pressure off the tuna population, which then starts growing. The cycle then repeats. — Modeling Life

The shark–tuna system is an example of a system with feedback. The tuna population has a positive effect on the shark population, while the shark population has a negative effect on the tuna population. This is an example of a negative feedback loop.

I plotted the change vector space for 10,000 cycles. The red circle indicates the starting point, and the green circle shows the state of the shark and tuna populations after 100 cycles. The vector field in the plot is moving in a counterclockwise direction. This is what’s causing the oscillations seen in the shark-tuna state space.

At the start, both the shark and tuna populations are decreasing. This is reflected in the vector field moving diagonally down and to the left. The vector field then shifts to moving left and up, indicating that the tuna population is increasing while the shark population is still declining.

Next, the vector field moves up and to the right, showing that both the shark and tuna populations are now increasing. Finally, the vector field shifts to moving to the right and down, signifying that the shark population is increasing while the tuna population is decreasing.

This cyclical pattern of the vector field demonstrates that the shark and tuna populations are oscillating over time, rather than following a linear trend. I highly recommend watching the shark tuna population animation.

Now, it’s time to explain why I multiplied the change vector (-5, -15) by a small value of 0.01 to get the changes of (-0.05, -0.15), and then used those scaled changes to update the shark and tuna state variables.

You can read my earlier post on How a Machine Learns to get a quick refresher on calculus. According to the Fundamental Theorem of Calculus (FTOC), given a function F, we can find its derivative F’ through differentiation, and through integration, we can recover the original function F from its derivative F’.

Consider the function 5X^3. The derivative of this function is 15X^2. The anti-derivative of 15X^2 is 5X^3+ C. Finding the derivative is easy, but integration is hard for real life examples. In fact, virtually none of the models we encounter in biology have a solution to the original function.

The change equations we derived tell us how the state variables S and T are changing with respect to time T. They’re S’(t) and T’(t). However, we don’t know the original functions S(t) and T(t). There’s no way to know them. So how were we able to plot the oscillating graph?

The answer lies in Euler’s method, a numerical technique discovered in 1768. By taking tiny steps, in our case using a time step of t = 0.01, we can approximate the original functions S(t) and T(t) even though we don’t know their exact equations.

Euler’s method consists in making ∆t very small but not zero. Starting at our initial point, we will follow its change vector for a very short time, specifically ∆t. We then find the change vector associated with the new point and follow it for the same very small time interval. Doing this over and over gives us a good approximation to the red curve, especially if we choose our ∆t to be very small. — Modeling Life

There’s a scene from the movie Hidden Figures where Taraji P. Henson played the role of Katherine Johnson, an American mathematician and NASA employee whose orbital mechanics calculations were critical to the success of the first and subsequent U.S. crewed spaceflights. In this scene, Katherine uses Euler’s method to calculate the paths for spacecraft to orbit the Earth and land on the Moon.

Back then there was an army of people called “computers” and 98% of them were women. It was simply not feasible to manually perform the thousands of computations required to calculate the spacecraft’s path. Katherine Johnson addressed the room full of women, saying, “Ladies, we are either computer programmers or we are unemployed.” They all learned to program, paving the way for the remarkable achievements that followed.

We mere mortals can’t reverse engineer all the complex functions that mother Nature has at her disposal to create all the living beings we see on Earth. Humans are clever. We don’t need to know the original function. All we need is to come up with a model by writing differential equations and use Euler’s method to figure out what the complex function is doing.

As Richard Feynman said, “You better learn Calculus, as it’s the language God speaks.” Modeling Life is one of the best introductory books for learning calculus. I will write a few more posts based on the insights I’ve gained from this book.