Population dynamics and functional programming

We don't need no iteration

Let’s play a game. We will write a model of population dynamics, and repeat it over a number of timesteps; then we will change the parameters, and draw a bifurcation diagram. These is the goal of the game. The only rule is as follows: if we write a loop, we lose. That’s right, we are going to do population dynamics using only functional programming. Because it is fun.

Setting up the model(s)

Let’s begin by defining two functions: one for exponential growth, and one for regulation. These are function generators, in that they will return a function, which can then be called with parameters.

R(r::Float64) = (n) -> n*r
K(k::Float64) = (n) -> n*(1.0-n/k)

We can create a growth model from these two components, by calling K first, then R. This is strictly useless, but a good opportunity to showcase the ∘ (\circ) operator, which in Julia works exactly like it works in the mathematical world, i.e. it composes functions.

growth = (R(2.3)∘K(1.0))

Let’s check what this does. First, it calls K to generate a function, which is $g(n) = 1-n\times K$, where $K$ is the carrying capacity, and $n$ is the population density. Then, it calls R to generate a function $f(n) = r\times Y$, where $Y$ is actually $g(n)$. So the final function is $f\circ g(n) = f(g(n))$, which is $r\times n\times(1 - n/K)$. This is the logistic model.

We can check the result of growth(1.0) (0.0), and growth(0.0) (0.0), as well as something above the carrying capacity, like growth(1.2) (-0.55).

This seems to be working! The next step is to actually run this a number of times, and to increment the population size. This can be achieved with the accumulate function. To accumulate something, we need an operator with two arguments, x[i] and x[i+1]. In our case, we want to add growth(n) to n, so we will not need the second argument, it’s simply here as a placeholder. Let us build the correct operator:

model = (current, _) -> current + growth(current)

Let’s apply this to a series of 50 elements, starting with 0.01. As the value of $r$ is large (2.3), this should lead to oscillations with a period of 2:

nt = accumulate(model, zeros(Float64, 50); init=0.01)
plot(nt, leg=false, c=:teal, frame=:origin)

Sweet - it works.

Bifurcation diagram

Now let’s ramp this up, and see if we can come up with a bifurcaton diagram without writing a for loop. To do this, we will need a sequence of values of $r$ to try, which is always the first step in a bifurcation diagram. As we have read our textbooks (or because we have used a very elegant approach to figure it out), we know that the system crashes above $r = 3.0$, so we will only go up until this point.

r_seq = 2.0:0.001:3.0

Now, we can map through these values, to create a series of models:

models = map(r -> (current, _) -> current + (R(r)∘K(1.0))(current), r_seq);

Every element in this array is a model, which can then be called within accumulate, as we did before. At this point, we can generate the bifurcation diagram:

# Generate the timeseries
raw_timeseries = map(m -> accumulate(m, zeros(Float64, 5000); init=0.01), models)

# Get the unique values across the latest 200 steps
endpoints = map(ts -> unique(ts[(end-200:end)]), raw_timeseries)

# Generate a list of the values of r for every simulation
rs = map(i -> fill(collect(r_seq)[i], length(endpoints[i])), eachindex(endpoints))

# Plot!
scatter(rs, endpoints, leg=false, c=:teal, msc=:transparent, msw=0.0, ms=1, alpha=0.2, grid=false, frame=:none)

This is exactly what we were expecting to see!

A conclusion?

This is, almost without a doubt, the most needlessly complicated way to solve this problem. But it goes to show that the tools of functional programming can be applied even for things that we would traditionally solve with a good old for loop. I particularly like the \circ notation, which is providing code that is ever closer to its mathematical meaning.