# 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.