Spruce Budworms

reproducing some results in May’77

Note

This is a demo combining the powers of Quarto with Entangled. Quarto is a utility for authoring notebooks in the Quarto variant of Markdown. Entangled is a tool that extracts code from codeblocks into runnable and reusable modules.

The codeblocks marked with >>> are evaluated through Jupyter by Quarto. Other code blocks either denote a named block with <<...>> quotation marks, or are written to an output file in the source tree by Entangled.

The result is a package that can be used by other projects without modification, so we’ve bridged the gap from a notebook to a package.

The source of this page is available at github.com/entangled/quarto-example.

This demo uses the numpy and matplotlib libraries.

>>>

import numpy as np
from matplotlib import pyplot as plt

%config InlineBackend.figure_format = 'svg'
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 14.0

The Model

We will reproduce some of the results found in the seminal paper by May (1977). This paper describes bifurcation in population models of the spruce budworm using this differential equation:

\[\frac{{\rm d} N}{{\rm d} t} = r N \left(1 - \frac{N}{K}\right) - \frac{cN^2}{H^2 + N^2},\]

where \(r\) is the intrinsic growth rate, \(K\) is the carrying capacity (related to the average leaf area), \(H\) is the characteristic budworm population (at which predation is saturated) and \(c\) is the level of predation. We store these parameters in a dataclass named BudwormModel:

<<budworm-model>>

@dataclass
class BudwormModel:
    growth_rate: float
    characteristic_population: float
    carrying_capacity: float
    predation: float | np.ndarray

    def set(self, **kwargs):
        return replace(self, **kwargs)

    <<budworm-methods>>

Here, the set method can be used to temporarily change any parameter value (by creating a copy). Let us first write down the differential equation in a form that Python understands, as a method of BudwormModel:

<<budworm-methods>>

def population_change(self, N):
    r = self.growth_rate
    K = self.carrying_capacity
    c = self.predation
    H = self.characteristic_population
    return r*N * (1 - N/K) - c*N**2 / (H**2 + N**2)

Let’s see how this function behaves for different values of predation. We may give the values of the other arguments in the form of a dictionary. We chose some values that should be reasonable.

>>>

from budworms.model import BudwormModel

model = BudwormModel(
    growth_rate=1.0,
    carrying_capacity=10,
    predation=2.0,
    characteristic_population=1.0)

The following plot shows the behaviour of this model for different levels of predation:

>>>

N = np.linspace(0, 8, 161)
fig = plt.subplot(111)  # obtain a figure object

for c in np.linspace(1.6, 2.8, 7):  # Loop over several values for `c`
    model.predation = c                        # change only the value of 'predation' in the model
    dN = model.population_change(N)            # apply our function for dN/dt
    fig.plot(N, dN, label='{:.03}'.format(c))  # plot

fig.axhline(0.0, c='k', linewidth=1)        # add a line on y == 0.0
fig.legend(loc='upper right', shadow=True)  # add a legend
fig.set_xlabel('N')                         # set labels
fig.set_ylabel('dN/dt')

plt.show()

Interesting things happen when derivatives vanish. Where lines cross the x-axis upward, we find unstable equilibria, and similarly, downward crossing gives a stable equilibrium. This plot tells us to expect interesting things (i.e. bifurcation) near \(c=1.8\) and \(c=2.6\). For those values the vanishing first derivative coinside with vanishing second derivative.

We can tell more about the situation at these special points by making a contour plot of \({\rm d}N/{\rm d}t\) with varying \(N\) and \(c\). To make matters easier we can get an analytic expression for the roots

\[\frac{dN}{dt} = 0,\]

which is solved by

\[N = 0\]

and

\[c = r\left(\frac{H^2}{N} - \frac{H^2}{K} + N - \frac{N^2}{K}\right).\]

<<budworm-methods>>

def critical_predation(self, N):
    H = self.characteristic_population
    K = self.carrying_capacity
    r = self.growth_rate
    return r*(H**2/N + N - H**2/K - N**2/K)

This expression implicitely tells us what the equilibrium values for \(N\) are, given \(c\). To find the tipping point values for \(c\), we should see where the derivatives vanish, giving

\[\frac{{\rm d} c}{{\rm d} N} = r\left(-\frac{2}{K} N^3 + N^2 - H^2\right) = 0.\]

This cubic polynomial can be solved using Cardano’s formula (for reference, see Wikipedia. We’ve enclosed an implementation of Cardano’s formula in Section 3.

<<budworm-methods>>

def critical_population(self):
    H = self.characteristic_population
    K = self.carrying_capacity
    return cubic_roots(-2/K, 1, 0, -H**2)

>>>

N_critical = np.sort(np.array(
    [e.real for e in model.critical_population()]))

print("critical N:        ", N_critical)
print("critical predation:", model.critical_predation(N_critical))

critical N:         [-0.919089   1.1378052  4.7812838]
critical predation: [-2.19159537  1.7872302   2.60436517]

The negative values are not real in the physical sense of the model, so we’ll ignore them. The values we find for \(c = 1.787\) and \(c = 2.604\) do match our expectations.

We’re now ready to combine all this information into a diagram, similar to Fig. 6 in May (1977).

>>>

fig = plt.subplot(111)

# add contour plot
N_cntr = np.linspace(-1, 10, 100)
c_cntr = np.linspace(1.6, 2.8, 100)
dN_cntr = model.set(predation=c_cntr[None,:]).population_change(N_cntr[:,None])
fig.contourf(c_cntr, N_cntr, dN_cntr, levels=np.linspace(-2, 2, 9), cmap='RdYlBu')

# emphasize the 0-contour and tipping points
N1 = np.linspace(0.01, N_critical[1], 50)
N2 = np.linspace(N_critical[1], N_critical[2], 50)
N3 = np.linspace(N_critical[2], 9.0, 50)
cn = model.critical_predation

fig.plot(cn(N1), N1, linewidth=2, c='k', label='stable equilibrium')
fig.plot(cn(N2), N2, linewidth=2, c='k', ls='dashed', label='unstable equilibrium')
fig.plot(cn(N3), N3, linewidth=2, c='k')
fig.plot(cn(N_critical), N_critical, 'o', ms=8, c='k')

# add quivers
quiver_mesh = np.meshgrid(
    np.linspace(1.65, 2.75, 20),
    np.linspace(-1, 8.5, 15))
N_q = quiver_mesh[1]
c_q = quiver_mesh[0]
quiver_dN = model.set(predation=c_q).population_change(N_q)
fig.quiver(c_q, N_q, 0, quiver_dN, width=0.0015, scale=20)

# decoration
fig.set_xlabel('predation (c)')
fig.set_ylabel('population (N)')
fig.set_ylim(0, 8.5)
fig.set_xlim(1.6, 2.8)
fig.legend(loc='upper right', shadow=True)

plt.show()

Packaging

In order to package our results we use the uv package manager. To create a Python package we need to create a __init__.py file in the package directory (here src/budworms).

src/budworms/__init__.py

# empty file

src/budworms/model.py

from dataclasses import dataclass, replace
from .cubic_roots import cubic_roots
import numpy as np

<<budworm-model>>

Cubic Roots

A cubic equation is solved using Cardano’s formula (for reference, see Wikipedia).

Given our equation,

\[ax^3 + bx^2 + cx + d = 0,\]

we first calculate,

\[\Delta_0 = b^2 - 3ac,\]

\[\Delta_1 = 2b^3 - 9abc + 27a^2d,\]

and,

\[C = \sqrt[3]{\frac{-\Delta_1 \pm \sqrt{\Delta_1^2 - 4\Delta_0^3}}{2}}.\]

Then we define the constant \(\zeta = -\frac{1}{2} + \frac{i}{2} \sqrt{3}\) (a rotation around 120 degrees). We can now find all three zeros as,

\[x = -\frac{1}{3a} \left(b + \zeta^k C + \frac{\Delta_0}{\zeta^k C}\right)\ \quad\textrm{for}\ k = 0, 1, 2.\]

src/budworms/cubic_roots.py

import numpy as np

def cubic_roots(a, b, c, d):
    """Compute the roots of the cubic polynomial :math:`ax^3 + bx^2 + cx + d`.

    :param a: cubic coefficient
    :param b: quadratic coefficient
    :param c: linear coefficient
    :param d: constant
    :return: list of three complex roots

    This function does not check if the found roots are real or complex.
    """
    if (a != 0): # Case: ax^3 + bx^2 + cx + d = 0
        delta_0 = b**2 - 3*a*c
        delta_1 = 2*b**3 - 9*a*b*c + 27*a**2*d
        C = ((delta_1 + np.sqrt(delta_1**2 - 4*delta_0**3 + 0j)) / 2)**(1/3)
        zeta = -1/2 + 1j/2 * np.sqrt(3)
        return [-1/(3*a) * (b + zeta**k * C + delta_0 / (zeta**k * C))
                for k in [0, 1, 2]]

    elif (b != 0): # Case: bx^2 + cx + d = 0
        delta = c**2 - 4*b*d
        return [(-c + k*np.sqrt(delta + 0j))/(2*b)
                for k in [-1, 1]]

    elif (c != 0): # Case: cx + d = 0
        return -d/c

    elif (d != 0): # Case: d != 0 (without solution)
        return np.nan

    else: # Case: d = 0 (with trivial solution)
        return 0

References

May, Robert M. 1977. “Thresholds and Breakpoints in Ecosystems with a Multiplicity of Stable States.” Nature 269: 471–77. https://doi.org/10.1038/269471a0.