Why is there a limit to growth in the logistic model?

A visual interpretation

Attention conservation notice: this entry is the aggregation of two drafts of lecture notes about using Julia for ecological modeling, each chapter being one lab session. Any feedback is particularly welcome, as the goal is to eventually turn these notes in a companion textbook. This entry is essentially two such modules stitched together, so it’s a little long. The point of these two modules is (i) to have a look at the bifurcation diagram of the logistic growth model, (ii) to see how we can build and read a cobweb plot, and (iii) how we can use the cobweb plot to develop an intuition about finding the maximal value of the growth rate. This uses base Julia, and a little bit of Symbolics to get the derivative for us.

using Plots
using Symbolics

The logistic map

In this section, we will have a look at (possibly) the most important single-population model from population biologist: logistic growth. In the logistic growth model, the dynamics of populaton growth are entirely governed by two parameters: its growth rate $r$, and its carrying capacity $K$. The models makes the assumption that all individuals have the same average number of offspring from one generation to the next, and that this number decreases when the population becomes more crowded, i.e. when $N(t) \rightarrow K$.

Defining the model

The expected number of offpsrings for each parent individual in the population is given by $1+r×(1-\frac{N(t)}{K})$; we can check that the population remains the same, i.e $N(t+1) = N(t)$ if the expected number of offpsrings is 1, which can happens when $r=0$, or when $N(t)=K$. In any case, every individual in the population will contribute the same number of offpsring, and so we can write the recursion equation of the logistic growth model as $N(t+1) = N(t)×(1+r×(1-\frac{N(t)}{K}))$. This function has two parameters ($r$, $K$), and one variable ($N(t)$).

From the above equation, we can calculate the expression of the change in population at every time step, $ΔN = N(t+1)-N(t) = N(t)×r×(1-N(t)/K)$. The population will grow when $ΔN > 0$, and shrink when $ΔN < 0$; as $r$ is a positive constant, and $N(t)$ is a positive (or null) population size, this means that the sign of $ΔN$, and therefore the growth of the population, is determined by the sign of $1-N(t)/K$.

timesteps = 100
N = zeros(Float64, timesteps + 1);

We need to initialize this array, so that N(t₀) has a value. The notation we use here refers to the first position as 1, but it is actually valid Julia notation to use N[begin].

N[1] = 0.1
0.1

One might argue that N is not a very explicit variable name; this is correct. Nevertheless, N is equivalent to $N$ in our mathematical notation, and when going back and forth between the model-as-mathematics and the model-as-computer-code, it helps decrease the cognitive load to have the same symbol mean the same thing. In short, we will sacrifice some programming good practices (“give variable explicit names”) in the name of consistency.

Simulating the model

Simulating the discrete-time logistic growth will be a simple question of iterating over values in the population size array, and calculating the change, then adding it to the next position. Because we know what happens at $t_0$, which is the position 1, we can either chose to iterate from 2 until the end, or from 1 until one before the end. As we expressed our model in terms of $N(t)$ and $N(t+1)$, rather than $N(t-1)$ and $N(t)$, we will start iterating at 1.

We do not, currently, have assigned values to either r or K. We will do this inside the loop, in order to not have parameter values in the global scope; after this illustration, we will write our model as a function, and give r and K default values.

for t in 1:timesteps
    local r, K = (0.2, 1.0)
    ΔN = N[t] * r * (1.0 - N[t] / K)
    N[t + 1] = N[t] + ΔN
end

We can now visualize the timeseries of this simulation – note that we specificy the limits on the time axis, because we start at 0, and end at the value determined by the timesteps variable.

scatter(0:timesteps, N; color = :black)
hline!([1.0], c=:grey, ls=:dash)
xaxis!("Time")
yaxis!("N(t)")

As we could expect, the population increases steadily, until it reaches a value given by $K$, here set at K = 1.0 and indicated by the horizontal dashed line. But what if we wanted to change $r$? Would the system behave in the same way? In order to make out work easier, which is to say more re-usable, we will write the model as a function:

function discrete_time_logistic_growth(
    N₀::Float64; tₙ::Integer=100, r::Float64=0.2, K::Float64=1.0
)
    N = zeros(Float64, tₙ + 1)
    N[1] = N₀
    for t in 1:tₙ
        ΔN = N[t] * r * (1.0 - N[t] / K)
        N[t + 1] = N[t] + ΔN
    end
    return (0:tₙ, N)
end
discrete_time_logistic_growth (generic function with 1 method)

Let’s take this function for a spin! We will use a specific syntax (x, y = z), in which the various elements of z are “unpacked” in x and y. We will also wrap everything within a begin ... end block, which is not strictly speaking necessary, but not a bad practice as it shows how we are treating the series of instructions as a single, cohesive unit. Note that we are using the keyword for the discrete_time_logistic_growth function to set the maximal timestep to 50, and the growth rate to a rather high 2.5.

begin
    t, N = discrete_time_logistic_growth(0.1; r=2.5, tₙ=50)
    scatter(t, N; color = :black)
    hline!([1.0], c=:grey, ls=:dash)
    xaxis!("Time")
    yaxis!("N(t)")
end

Compared to the previous figure, the model shows a qualitatively different behavior; the behavior of two models is said to differ qualitatively when, in order to describe what you see, you cannot simply rely on words like “lower” or “earlier” or “more frequently”; in this case, we could for example say that whereas with $r=0.1$ the population reached a value that looked like it was $K$, using $r = 2.5$ leads to oscillations around the value of $K$.

The bifurcation diagram

From the previous section, it looks like the behavior of the model depends on the value of $r$ – in the next chapter, we will explore a method to identify the maximal value of $r$ after which the population is guaranteed to go extinct, but for now we will accept that it is $r = 3$. One question we can ask is: what are the values taken by the population size after a sufficiently long period of time has passed? The purpose of letting enough time pass is to make sure that, for slow growth dynamics (i.e. low $r$), we give the model enough time to reach the equilibrium.

The figure we will use to visualize this is called a bifurcation diagram; building it is relatively simple: we take all values of $r$, get the last 1000 (or any other arbitrary constant) values of $N(t)$, and plot the relationship between the two.

all_r = Float64[]
all_N = Float64[]

for r in LinRange(1.8, 3.0, 500)
    t, N = discrete_time_logistic_growth(0.1; tₙ=5_000, r=r)
    last_N = unique(N[(end - 1000):end])
    last_r = fill(r, length(last_N))
    append!(all_N, last_N)
    append!(all_r, last_r)
end

scatter(all_r, all_N; color=:black, dpi=600, ms=1, alpha=0.1, lab="")
yaxis!("N(t)")
xaxis!("r")

Revisiting the logistic map

When models are simple enough, here meaning that they do not have a large number of variables, it is possible to develop often deep intuitions about their behavior simply by visualizing some of their properties. In this chapter, we will look at discrete and continuous time models, and explore some of Julia’s automatic differentiation capacities to figure out when derivatives change.

In the previous section, we spent a lot of time exploring the behavior of the logistic map. We will now see how much of this behavior we could have expected form very simple visualisations. To begin with, we will write a function that returns a function. Why? There are two reasons to this choice. The first is very pragmatic: it will make some of the following code easier to read. The second is to demonstrate some functionalities of Julia: functions are “first-class citizens”, and can be created and returned just like any other object.

function logistic_growth_generator(r::T, K::T) where {T<:Float64}
    return (N) -> N * (1 + r * (1 - N / K))
end
logistic_growth_generator (generic function with 1 method)

How do we use this function? Calling logistic_growth_generator will return a “generic” function, which accepts one argument (N; the population size), but where r and K are built-in. We can use this to ask questions like, what will be the population size at the next generation for a population of $N = 0.8$, with $r = 1.2$ and $K = 0.9$?

logistic_growth_generator(1.2, 0.9)(0.8)
0.9066666666666667

We can also use this function to ask the same question of a range of values:

logistic_growth_generator(1.2, 0.9).(0.0:0.3:0.9)
4-element Vector{Float64}:
 0.0
 0.54
 0.84
 0.9

This is all we need to start building a plot. The specific type of plot we usually build for a single-variable discrete-time model is a “cobweb plot”, which represents the space of $N(t)$ and $N(t+1)$; the two elements are the 1:1 line, and the line defined by the recursion equation. Whenever the recursion equation intersects the 1:1 line, we know that $N(t) = N(t+1)$, or in other words $\Delta N = 0$, i.e. the system is at an equilibrium.

The cobweb plot

Although we may want to get a collection of values, we can use the fact that Julia is able to natively plot functions. The idea of writing a single-variable function will therefore make this code very simple to write, especially if we want to add multiple functions!

plot((x) -> x; xlim=(0, 1.6), ylim=(0, 1.6), aspectratio=1, lab="", c=:grey, ls=:dash, frame=:box, dpi=600)
plot!(logistic_growth_generator(0.5, 1.0); lab="r = 0.5", c=:black)
plot!(logistic_growth_generator(1.5, 1.0); lab="r = 1.5", c=:teal)
plot!(logistic_growth_generator(2.5, 1.0); lab="r = 2.5", c=:orange)
xaxis!("N(t)")
yaxis!("N(t+1)")

Let us focus on an example, and try to follow the trajectory of a population along this plot – the origin of the name “cobweb” should become apparent. We will pick $r = 2.35$ and $K = 1.0$, with $N(t_0) = 0.04$. How does this population move in this diagram? We start at the point $(N(t), N(t))$, then move up to $(N(t), N(t+1))$, then project on the diagonal to $(N(t+1), N(t+1))$, and then finally we start the cycle again.

begin
    𝑓 = logistic_growth_generator(2.35, 1.0)
    N = 0.04
    plot(
        (x) -> x;
        xlim=(0, 1.6),
        ylim=(0, 1.6),
        aspectratio=1,
        lab="",
        c=:grey,
        ls=:dash,
        frame=:box,
        dpi=600,
    )
    plot!(𝑓; lab="", c=:black, lw=2)
    for t in 1:200
        global N
        N_next = 𝑓(N)
        plot!([N, N, N_next], [N, N_next, N_next]; lab="", c=:darkgrey, alpha=0.6)
        N = N_next
    end
    xaxis!("N(t)")
    yaxis!("N(t+1)")
end

This diagram is an example of a “cobweb” plot, named for the fact that the grey line, representing the jumps made by the system across discrete time, look like a cobweb. Each jump obeys the same sequence: move to the line of the recursion equation, then horizontally to the 1:1 line, then vertically to the line defined by the recursion equation.

We can see the system initial creeping up towards the equilibrium point, but the high growth rate means that it will get overshot, and ultimately the system reaches a series of orbits around this point. Is there a deep biological intuition hidden in this figure? None that we could not get from the bifurcation diagram. But there is an intutition that can lead us to a mathematical result, specifically the maximal value of $r$ for which the system will not go extinct!

Figuring out the maximal growth rate

Notice one thing on this plot: when the system reaches the maximal value of the recursion equations, i.e. the top of the parabola, it gets projected on the 1:1 line, at a point that can be above further right (on the $N(t)$ axis) than the equilibrium. It then gets projected down towards the curve for the recursion equation once more; if the value on the $N(t+1)$ axis of this projection is negative, then the population goes extinct. Can we figure out this threshold?

Before we do this, let’s think some more about the point of inflection of the curve corresponding to the recursion equation. For high values of $r$, it does happen to the left of the equilibrium, but this is not true for loe values of $r$. In order to figure out where this maximum is, we can for example calculate $N(t+1)$ for a lot of values of $N(t)$, and pick the maximum:

𝑓 = logistic_growth_generator(2.35, 1.0)
Nt = 0.0:0.001:1.6
Nt[Nt .|> 𝑓 |> argmax]
0.713

This gives the value of $N(t)$ for which $N(t+1)$ is the largest, and of course we can get the value of this largest possible $N(t+1)$ by

𝑓(Nt[Nt .|> 𝑓 |> argmax])
1.19388285

Take a moment to read this on the figure, and verify that it matches the visual representation of the system. When this is done, we can verify that $N(t+2)$ is still positive (the system for the parameters we have chosen does not go extinct), using

𝑓(𝑓(Nt[Nt .|> 𝑓 |> argmax]))
0.649920337618312

Note that the notation $g(g(x))$ is also $(g∘g)(x)$, and Julia does allow this exact notation:

(𝑓∘𝑓)(Nt[Nt .|> 𝑓 |> argmax])
0.649920337618312

With this in hands, we can start thinking about the value of $r$ for which we will start a shift from a steady approach of the equilibrium to oscillations (that may or may not be sustained): it hapens when the maximum of $𝑓$ is larger than $K$. We can check that $f(x)$ has an extremum at $a$ if $f'(a)=0$ and $f'(x)$ is defined at $a$. We could easily calculate this derivative ourselves, but we can also rely on Julia’s Symbolics package to do it for us – we will prefix the variables with s to mark that these represent symbols, so sK means “the symbol $K$”.

@variables sr sK sN
sNt = sN*(1+sr*(1-sN/sK))
D = Differential(sN)
∂Nt = simplify(expand_derivatives(D(sNt)))
(sK + sK*sr - 2.0sN*sr) / sK

Since $K > 0$, the maximum of the recursion equation is reached when $K(1+r)=2rN(t)$, which is $N(t) = K(1+r)/(2r)$. We are now much closer to figuring out the maximal value of $r$! We know which value of $N(t)$ gives the maximum $N(t+1)$, and so we need to figure out when this $N(t+1)$, when projected on the 1:1 line, and then on the recursion equation curve once again, will be negative. But the last step is actually contained in the first one – looking back at the cobweb plot, all that is needed for $N(t+2)$ to be negative is for $N(t+1)$ to be larger than the value of $N(t)$ for which the recursion equation is equal to 0. This last part is easily solved, discarding the obvious $N(t) = 0$, and gives a solution of $N(t) = K(1+r)r^{-1}$. Try modifying the cobweb plot to visualize the vertical lines corresponding to $N(t)$ giving the maximal $N(t+1)$, and $N(t+1)=0$; this is a good way to confirm that your calculations are correct.

And just like this, our problem is almost solved! We need to calculate the value of $N(t+1)$ when $N(t)$ gives the maximal value:

Nmax = sK*(1+sr)/(2*sr)
Ntmax = simplify(substitute(sNt, Dict(sN => Nmax)))
(sK*((1//2) + (1//2)*sr)*(1 + sr)) / (2sr)

For this value (let’s call it $x$) to be larger than $y = K(1+r)(1/r)$, we can check the conditions under which $x/y > 1$. This turns out to simplify quite nicely:

Ntmax / (sK*(1+sr)*(1/sr))
((1//2) + (1//2)*sr) / 2

After a bit of re-organisation (note how $K$ disappears!), we are left with the solution to our problem: $r > 3$. In short, for any value of $r$ such that $0 < r ≤ 3$, the discrete-time logistic growth model will persist. It is important to note that we reached this conclusion by looking at the cobweb plot: there is a point past which the projection is going to bring us into negative values of $N(t)$. Since this is impossible, we thought about the positions of various points on this diagram that would get us in this situation, and this ended up having a simple solution, namely that $r$ must be at most 3.