Well, I feel like I should apologize for such a long delay between posts. It’s been a crazy summer that has included some vacation time plus an overseas trip to a conference. Regardless, I’m finally back in the swing of things and ready to write up Part 3! To recap: Lesson 1 and Lesson 2 looked at how ODEs are solved numerically and how higher order solutions are more accurate than naive implementations. Today we’ll look at two simulations of living systems (Lotka-Volterra and SIR). Finally, the series will conclude with a post on model fitting and a post about chaotic systems.
Previous posts explained how numerical solutions work and how Matlab will perform the calculations for you automatically. Today, we’ll put that knowledge to good use. We’ll start with a simple Lotka-Volterra predator/prey two-body simulation. There are two populations in question: the predators and the prey. Prey multiply exponentially, similar to our exponential example in the previous lessons. The difference is that prey are also killed off by the predators at a rate directly proportional to both the predator and prey population. Predators are dependent on prey for sustenance and thus grow at a rate dependent on both the predator and prey population. Predators also die off over time due to age.
This gives us two populations each with two separate terms to add to our equation. Because both predator and prey start with ‘p’, we will use the variable x to refer to the prey population and the variable y to refer to the number of predators. Given that both x and y are functions of time, we have the following system of ODEs:
Let’s break things down. The top equation defines how the population of prey, x, changes in relation to the predator and prey populations. The first term, αx, is positive and therefore defines growth. The growth rate is proportional to x, so we know it is positive exponential growth (as seen in previous lessons). The α term is a constant that controls how fast the exponential growth occurs. The second term, -βxy, is negative and thus describes death. The rate is proportional to both x and y. This is because more prey makes it easier for predators to find prey and also because an increase in predators makes it harder for the prey to hide. In contrast, we assume an increase in predators does not affect the birth rate of the prey and thus there is no y in the growth term. Thus, βxy defines the rate at which predators find and kill prey, at an efficiency controlled by the constant β.
The second equation defines the change in the predator population. The first term, δxy is positive and defines growth. Predators rely on prey for sustenance and thus grow proportionally to both their own population and the prey population at a rate controlled by δ. For the second term, because prey can’t kill predators, the predator population only decays proportionally to their own numbers at a rate determined by γ.
It is important to note that both equations have an xy term. This means that both equations are non-linear and thus very difficult to deal with analytically. This is one of the reasons we will turn to numerical integration.
Now, let’s look at the expected dynamics before we actually plug in the equation. We have two populations, so let’s look at the different combinations of high and low populations of each. We want to pick coefficients such that the model is stable. If both populations are low, there needs to be positive growth to restore balance. This growth should come from the prey population being able to multiply without fear of many predators. Thus, the α term must out-pace the β term for small values of x and y. If both values are high, there needs to be death to restore equilibrium. Thus the prey needs to die off due to a large predator population and the β term will dominate. If there is a lot of prey with few predators, the predator population should grow quickly. Conversely, if there are many predators and few prey, the predators population should diminish. If all of this happens, our system should be able to stay in some sort of periodic equilibrium. To this end, I will somewhat arbitrarily chose the constants to be: α = 1, β = .05, δ = .02, and γ = .5.
First, we will code up the equations as a matlab function. The difference between this code and the code from Lesson 2 is that now we are dealing with two variables. Thus, our x parameter is now a vector containing both the prey and predator populations.
% lotka_volterra.m % % Imlements the Lotka-Volterra function % dx/dt = alpha x - beta xy % dy/dt = delta xy - gamma y % % Inputs: % t - Time variable: not used here because our equation % is independent of time, or 'autonomous'. % x - Independent variable: this contains both populations (x and y) % Output: % dx - First derivative: the rate of change of the populations function dx = lotka_volterra(t, x) dx = [0; 0]; alpha = 1; beta = .05; delta = .02; gamma = .5; dx(1) = alpha * x(1) - beta * x(1) * x(2); dx(2) = delta * x(1) * x(2) - gamma * x(2);
Unless you plan on modifying your coefficients through external code it’s better to define them inside the function. Also note that x(1) refers to the prey and x(2) refers to predators.
Next, we need to use odeset to set our options. Note that now both populations 1 and 2 are non-negative.
% Set our preferences for ode45 % The default relative tolerance is 1e-3. % To set our output to non-negative, we provide an array listing % each population that we want to constrain. Since this example % has two populations, we pass the array [1 2] options = odeset('RelTol', 1e-4, 'NonNegative', [1 2]);
Finally, we use the ode45 function and plot out the results. We will initialize each population to 10 and will integrate from time 0 to 20. The time units here are completely arbitrary since we are using unitless coefficients.
% Use ode45 to solve our ODE % Place the time points in a vector 't' % Place the solution in a vector 'x' [t,x] = ode45('lotka_volterra', [0 20], [10 10], options); plot(t,x); legend('prey', 'predators');
Which gives us a final result pictured here:
Notice the nice cyclical behavior. At first, prey multiply rapidly and predators decline a bit. As the prey population soars, predators find more food and begin to grow as well. Eventually the predators eat too many prey and lose numbers to starvation. This lets the prey rebound and the cycle begins again.
Let’s look at one more example. An SIR model is an epidemiological example of an infection invading a population. In the beginning most people are healthy and the infection spreads slowly. As more people get sick, the infection begins to grow noticeably. Eventually, the population begins to recover and gains an immunity to the infection. Thus, the susceptible population disappears and is replaced by a resistant population. The three populations modeled in the standard SIR equation are susceptible, infected, and recovered.
The first thing to note is that there are only two terms. The first, βSI describes the rate of infection. This term scales positively with both the susceptible population and the infected population. The more people susceptible, the easier it is to find a new host. The more infected people, the more there are to spread the disease. Also note that the exact same quantity that is removed from the susceptible population is added to the infected population. This represents a direct state transition from one population to the other. Similarly, as infected people gradually recover at a rate δI, they are removed from the infected population and added to the recovered population. Finally, note that if you add all three equations together you get a net effect of (dS+dI+dR)/dt = 0. This shows that the total population is stable: there is no addition or removal of people, only changes in their state.
Let’s code this up! First, the function definition (I will set the constants to be β = .003 and δ = 1:
% sir.m % % Imlements a SIR infection model % dS/dt = -beta SI % dI/dt = beta SI - delta I % dR/dt = delta I % % Inputs: % t - Time variable: not used here because our equation % is independent of time, or 'autonomous'. % x - Independent variable: this contains three % populations (S, I, and R) % Output: % dx - First derivative: the rate of change of the populations function dx = sir(t, x) dx = [0; 0; 0]; beta = .003; delta = 1; dx(1) = -beta * x(1) * x(2); dx(2) = beta * x(1) * x(2) - delta * x(2); dx(3) = delta * x(2);
We will use the same options, but add one extra non-negative population (as we now have three).
options = odeset('RelTol', 1e-4, 'NonNegative', [1 2 3]);
Finally, we will model an infection of a population of 1000 healthy individuals and 1 sick person. The time interval will go from 0 to a unitless 10.
% Use ode45 to solve our ODE % Place the time points in a vector 't' % Place the solution in a vector 'x' [t,x] = ode45('sir', [0 10], [1000 1 0], options); plot(t,x); legend('S', 'I', 'R');
And here is our result:
Notice how the infection initially picks up steam and spreads quickly through the population but soon runs out of gas as the susceptible population quickly drops. The simulation will continue indefinitely as the infected population will slowly decay and the recovered population will asymptotically approach 1001.
I will end with one final point about running continuous simulations of discrete populations. ODE models are a very nice way to simulate complicated population dynamics. Unfortunately, they can only be at best a good approximation. Populations are not continuous. Through ODE modeling, it is possible to represent positive populations of sizes less than one. Our SIR model would work just as well if we started it with an infected population of one tenth a person. Our Lotka-Volterra system would still be cyclical even if our predator population dropped down to one one-hundredth of an individual. Clearly this is not realistic. The lesson here is that you must always take the time to evaluate your model from a rational perspective. Don’t fall into the trap of only looking at the quantitative dynamics. Before you use