Mathsy Coding

Logistic Growth

Working through the logistic growth model
Edits:
  • : Changing image and file paths
  • : Changing picture and file paths

Introduction and Explanation

Note: most of this is covered in Chapter 1 of (Murray, 2002).

The simplest model for the growth of a population of organisms is that given by linear growth. This model assumes the rate at which organisms grow / reproduce is proportional to the number that currently exist. For example, imagine a bacterial culture where each bacterium buds into two daughter. This can be represented as a differential equation.

dNdt=rN, N(0)=N0\frac{dN}{dt} = rN,\, N(0)=N_0

Where N(t)N(t) is the size of the population at a time tt, rr is the growth rate, and N0N_0 is the initial population size. The solution to this is given by N(t)=N0ertN(t) = N_0e^{rt}. I’ve written about this model before.

While this model has the virtue of simplicity, it is not particularly realistic. The size of the population, assuming that the growth rate rr is positive, quickly reaches ludicrous levels. As such, we need a better model!

One such model, which accounts for the depletion of resources, accumulation of wastes, and other such inhibitors of growth that appear as the population grows, is the logistic model. The reasoning behind this model is that in the absence of any inhibiting factors, the population will grow as per the linear model. However, as the population grows, there is a factor that accounts for the slowing growth caused by depletion of resources, &c. As a differential equation, this looks as follows:

dNdt=rN⏟Linear growth(1−NK)⏞Inhibiting factor\frac{dN}{dt} = \underbrace{rN}_{\text{Linear growth}}\overbrace{\left( 1 - \frac{N}{K} \right)}^{\text{Inhibiting factor}}

where rr is the growth rate (as in the linear model) and KK is the carrying capacity.

Stationary Points

From the above equation, the system has two stationary points: N=0N = 0 and N=KN = K. Intuitively, these make sense: if there is no initial population, we don’t expect spontaneous generation to create one. Similarly, if the population is exactly at the carrying capacity, then we expect it to remain constant. We can also see that if NN is slightly larger than KK, the derivative is negative - the population will decrease towards the carrying capacity. If NN is slightly smaller than KK, then the derivative will be positive, and the population will increase towards the carrying capacity.

Let’s take a look at the slope field for this system, settings r=1r=1 and K=100K=100:

library(ggplot2)

len <- 5
diff <- function(N) N * (1 - N / 100)

ts = seq(0, 10, by=1)
Ns = seq(0, 120, by=10)

df <- expand.grid(t=ts, N=Ns)
df$m <- diff(df$N)
df$theta <- atan(df$m)
df$xstart <- df$t - 0.5 * len * cos(df$theta)
df$ystart <- df$N - 0.5 * len * sin(df$theta)
df$xend <- df$t + 0.5 * len * cos(df$theta)
df$yend <- df$N + 0.5 * len * sin(df$theta)

ggplot(df) +
    geom_segment(aes(x=xstart, y=ystart, xend=xend, yend=yend)) +
    xlab("t") +
    ylab("N") +
    ggtitle("Slope Field for Logistic Growth") +
    theme(plot.title=element_text(hjust=0.5))

Slope field for logistic growth

This makes it look like the stationary point N=0N=0 is unstable and N=KN=K is stable. Now let’s justify that algebraically. To do that, we’ll linearize the equation by taking a Taylor series about each of the stationary points and dropping all but the first two terms. When we do this, we find that a system (in general, not just this one) is stable if the second derivative at that point is negative and unstable if it is positive.

Say that N∗N^{*} is a stationary point of dNdt=f(N)\frac{dN}{dt} = f(N) so that f(N∗)=0f(N^{*}) = 0. Let’s consider δ(t)=N(t)−N∗\delta(t) = N(t) - N^{*}, the difference between the current population and the stationary point. Then

dδdt=dNdt−dN∗dt⏟0, since N∗ is constant=dNdt=f(N)=f(δ+N∗)\begin{align*} \frac{d\delta}{dt} &= \frac{dN}{dt} - \underbrace{\frac{dN^{*}}{dt}}_{\text{0, since $N^{*}$ is constant}} \\ &= \frac{dN}{dt} \\ &= f(N) \\ &= f(\delta + N^{*}) \\ \end{align*}

Then we can expand the Taylor series for dδdN\frac{d\delta}{dN} around N∗N^{*}

dδdt=f(N)=f(N∗)+f′(N∗)(N−N∗)+f′′(N∗)(N−N∗)2+…≈f(N∗)+f′(N∗)(N−N∗)Linear approximation=0+f′(N∗)(N−N∗)N∗ is a stationary point=f′(N∗)(δ+N∗−N∗)=f′(N∗)∗δ\begin{align*} \frac{d\delta}{dt} &= f(N) \\ &= f(N^{*}) + f^{\prime}(N^{*})(N - N^{*}) + f^{\prime\prime}(N^{*})(N - N^{*})^2 + \dots \\ &\approx f(N^{*}) + f^{\prime}(N^{*})(N - N^{*}) & \text{Linear approximation} \\ &= 0 + f^{\prime}(N^{*})(N - N^{*}) & \text{$N^*$ is a stationary point} \\ &= f^{\prime}(N^{*})(\delta + N^* - N^{*}) \\ &= f^\prime(N^*)\ast \delta \end{align*}

So then we have that, at least in a small region around the equilibrium point, dδdt=f′(N∗)∗δ→δ(t)=Aef′(N∗)t\frac{d\delta}{dt} = f^\prime(N^*)\ast \delta \to \delta(t) = Ae^{f^\prime(N^*)t}. From this, it is clear that if f′(N∗)>0f^\prime (N^*) > 0, then this is exponential growth - the difference between the population and the equilibrium population will grow rapidly, and so the solution is unstable. On the other hand, if f′(N∗)<0f^\prime (N^*) < 0, then this is exponential decay, and so the difference between the population and the equilibrium population will go to zero - the solution is stable.

For our particular system, f′(N)=ddtrN(1−NK)=r−2rKNf^\prime(N) = \frac{d}{dt} rN \left( 1 - \frac{N}{K} \right) = r - \frac{2r}{K}N. So, for the two equilibrium points:

Thus, we’ve confirmed algebraically what our intuition and the slope field diagram led us to expect.

Solution

Luckily for us, we can solve this exactly by using a separation of variables.

dNdt=rN(1−NK)dNN(1−NK)=rdt∫dNN(1−NK)⏟LHS=∫rdt\begin{align*} \frac{dN}{dt} &= rN\left( 1 - \frac{N}{K} \right) \\ \frac{dN}{N\left( 1 - \frac{N}{K} \right)} &= r dt \\ \underbrace{\int \frac{dN}{N\left( 1 - \frac{N}{K} \right)}}_{\text{LHS}} &= \int r dt \tag{1} \end{align*}

To solve the left hand side (LHS), we can find the partial fraction decomposition of the integrand:

1N(1−NK)=AN+B1−NK1=A(1−NK)+BN1=A+N(B−AK)∴A=1, B=1K\begin{align*} \frac{1}{N\left( 1 - \frac{N}{K} \right)} &= \frac{A}{N} + \frac{B}{1-\frac{N}{K}} \\ 1 &= A\left(1 - \frac{N}{K}\right) + BN \\ 1 &= A + N\left(B - \frac{A}{K} \right) \\ \therefore A=1,\ &B=\frac{1}{K} \end{align*}

So now, continuing on from (1),

∫dNN(1−NK)=∫rdt∫1N+1/K1−NKdN=∫rdt∫1N+1K−NdN=∫rdtln⁡∣N∣−ln⁡∣K−N∣=rt+Cln⁡∣NK−N∣=rt+CNK−N={AertN<K−AertN>K\begin{align*} \int \frac{dN}{N\left( 1 - \frac{N}{K} \right)} &= \int r dt \\ \int \frac{1}{N} + \frac{1/K}{1 - \frac{N}{K}}dN &= \int r dt \\ \int \frac{1}{N} + \frac{1}{K - N}dN &= \int r dt \\ \ln |N| - \ln |K - N| &= rt + C \\ \ln \left| \frac{N}{K - N} \right| &= rt + C \\ \frac{N}{K - N} &= \begin{cases} Ae^{rt} & N < K \\ -Ae^{rt} & N > K \\ \end{cases} \\ \end{align*}

For the case where N<KN < K:

N(1+Aert)=KAertN=KAert1+AErt=Kert1A+ert\begin{align*} N(1 + Ae^{rt}) &= KAe^{rt} \\ N &= \frac{KAe^{rt}}{1 + AE^{rt}} \\ &= \frac{Ke^{rt}}{\frac{1}{A} + e^{rt}} \end{align*}

By the same process, when N>KN > K, N(t)=Kert−1A+ertN(t) = \frac{Ke^{rt}}{-\frac{1}{A} + e^{rt}}.

If we have that N(0)=N0N(0) = N_0, then for the case where N<KN < K,

N0=N(0)=K1A+11A+1=KN01A=KN0−1\begin{align*} N_0 &= N(0) = \frac{K}{\frac{1}{A} + 1} \\ \frac{1}{A} + 1 &= \frac{K}{N_0} \\ \frac{1}{A} &= \frac{K}{N_0} - 1 \\ \end{align*}

And so putting this all together, we have that

N(t)=KertKN0−1+ert=N0KertK−N0+N0ert=N0KertK+N0(ert−1)\begin{align*} N(t) &= \frac{Ke^{rt}}{\frac{K}{N_0} - 1 + e^{rt}} \\ &= \frac{N_0 K e^{rt}}{K - N_ 0 + N_0 e^{rt}} \\ &= \frac{N_0 K e^{rt}}{K + N_0 ( e^{rt} - 1)} \end{align*}

Similar logic for the case N>KN > K gives the identical solution (which feels slightly miraculous).

Nondimensional Form

The last thing I’d like to do is to write down the nondimensional form of the equation. To do that, we’ll first substitute in values of the variables in our system (NN and tt) with normalized versions, N∗N^* and t∗t^*, each of which is a multiple of the original variables. If we choose the multiplying constants correctly, we can remove all constants in the equation (like rr and KK), leaving us with something just in terms of the variables themselves.

To see how this works, let’s start with a simple example: linear growth.

dNdt=rN\frac{dN}{dt} = rN

We’ll make the substitutions N=AN∗N = AN^*, t=bt∗t = bt^*. Then

dNdt=rNd(AN∗)d(Bt∗)=r(AN∗)ABdN∗dt∗=rAN∗dN∗dt∗=rB⏟Should =1N∗∴rB=1→B=1r\begin{align*} \frac{dN}{dt} &= rN \\ \frac{d(AN*)}{d(Bt^*)} &= r(AN^*) \\ \frac{A}{B} \frac{dN^*}{dt^*} &= rAN^* \\ \frac{dN^*}{dt^*} &= \underbrace{rB}_{\text{Should $=1$}}N^* \\ \therefore rB &= 1 \to B = \frac{1}{r} \\ \end{align*}

Thus, with the variable substitutions N∗=NN^* = N and t∗=rtt^* = rt, we have our nondimensional form. To verify:

dN∗dt∗=1rdNdt=1rrN=N=N∗∴dN∗dt∗=N∗\begin{align*} \frac{dN^*}{dt^*} &= \frac{1}{r}\frac{dN}{dt} \\ &= \frac{1}{r} rN \\ &= N \\ &= N^* \\ \therefore \frac{dN^*}{dt^*} &= N^* \\ \end{align*}

In this case, we can think of our new variable t∗t^* are representing the ‘characteristic time’ of the system - that is, we measure time in terms of the time it takes to grow by a factor ee. In some way, this is the natural way to measure time in this system.

Now let’s try the same process using our actual equation. Again, we’ll try the substitutions N=AN∗N = AN^*, t=Bt∗t = Bt^*. Then

dNdt=rN(1−NK)ABdN∗dt∗=r(AN∗)(1−AN∗K)dN∗dt∗=rB⏟=1N∗−rABK⏟=1(N∗)2\begin{align*} \frac{dN}{dt} &= rN\left( 1 - \frac{N}{K} \right) \\ \frac{A}{B} \frac{dN^*}{dt^*} &= r (AN^*)\left( 1 - \frac{AN^*}{K} \right) \\ \frac{dN^*}{dt^*} &= \underbrace{rB}_{=1}N^* - \underbrace{\frac{rAB}{K}}_{=1}(N^*)^2 \\ \end{align*}

So B=1rB = \frac{1}{r} and A=KA = K, giving us the subsitutions N∗=NKN^* = \frac{N}{K} and t∗=rtt^* = rt. When we make this subsitution, we have the nondimensional form

dN∗dt∗=N∗(1−N∗)\frac{dN^*}{dt^*} = N^*(1 - N^*)

Again, we can think of these as having physical significance. These subsititions mean that we are measuring time in terms of the characteristic growth rate time, and we are measuring population in terms of the carrying capacity.

Conclusion

The logistic growth model, and the curve that arises from it naturally, are fascinating to look at. In this post, I’ve just scratched the surface - there are many other interesting things to learn about. The logistic function, as well as arising in various forms in many ecological systems, also arises in other, seemingly unrelated fields. For instance, it is often used in machine learning as well as in the ranking of individuals. Hopefully this helped to explain some of the properties of this system and the resulting function!

Sources and Futher Reading