Modeling Population Recovery Following an Environmental Disturbance

By Azmy S. Ackleh and Amy Veprauskas

What do populations of invasive and endangered species have in common? To find out, consider the following two contrasting scenarios.

Small mammals, such as kangaroo rats, are considered to be keystone species in many grassland and shrubland communities [7]. This means that their presence and densities help shape community composition. Because they are granivores, kangaroo rats significantly impact annual plants that serve as a resource base [4]. However, these communities are subject to various disturbances—including habitat fragmentation, fires, and livestock grazing—that degrade habitat quality and regularly threaten keystone species. For kangaroo rats, dense covers of herbaceous nonnative plants magnify the effects of these disruptions [5]. This type of vegetation cover has been shown to affect population recovery following a disturbance.

In contrast, American bullfrogs are an invasive species that damage native fauna in habitats around the world [8]. The bullfrog tadpoles’ voracious appetites may dramatically reduce algae biomass—in turn reducing primary production and nutrient cycling—while adults compete with native species of birds, reptiles, amphibians, and fish for food sources [9]. As such, bullfrogs may have profound effects on native habitats, changing ecosystem structure and even causing local extinctions among native species. Control methods for American bullfrogs typically focus on the removal of tadpoles or adults from the population.

What similarities exist between these two scenarios? In addition to both populations substantially impacting overall community structure, these situations have two common components. Both cases involve the idea of a disturbance — such as habitat fragmentation or fires in the first case, or the intentional removal of individuals in the latter. There is also the concept of recovery. Recovery for an invasive species means that management strategies must be reapplied or modified. Recovery of an endangered species is the end goal.

Biological populations continually experience natural and anthropogenic disturbances—like hurricanes, fires, and chemical and noise pollution—that negatively influence their growth. From a management perspective, it is important to be able to quantify the way in which disturbances may affect a population’s dynamics over time. This knowledge can help set harvest or land use regulations, identify effective conservation approaches, or aid in the establishment of control measures for pest species. However, one must exercise caution when applying management strategies; in some cases, they may have unintended effects. For instance, studies have shown that the aforementioned strategy for bullfrog regulation might be an ineffective means of population control, and in some cases actually result in increased population sizes. Instead, removal of metamorphs in the fall may be more effective [8]. And while livestock grazing is generally believed to negatively impact some species’ environments, research has indicated that it can actually promote kangaroo rat recovery if it reduces vegetation coverage [6].

Mathematical modeling can serve as a complementary tool to experimental studies for understanding the implications of management or control strategies. It is both inexpensive and able to provide real-time management methods that do not require extended periods of data collection. However, model predictions rely heavily on the data’s accuracy and the assumptions used in model development. Nevertheless, models can provide useful insights—even when limited data is available—and help generate hypotheses that inform experimental designs. Here we present a general modeling approach for the study of population recovery. This approach is adaptable to various situations and may assist in the identification of effective control strategies.

How Can We Model Recovery?

Depending on the population under consideration, a population’s recovery may take many forms. If we wish to mathematically model recovery, we must first define the way in which a population’s size changes over time. To do so, we describe a population using a matrix model that allows us to distinguish between individuals in different developmental stages. 

Consider a female population that is divided into \(m\) stages. Denote the densities of these stages at time \(t\) with \(\textbf{n}(t):=[n_1(t),n_2(t),\ldots,n_m(t)]^{\intercal}\), where \(\intercal\) signifies the transpose of a vector. Let \({\textbf{A}}[\boldsymbol{\theta}({\epsilon}(t),\textbf{n}(t)),\epsilon(t)]\) be the projection matrix for the population at time \(t\) that describes individuals’ transitions between the different stages. This matrix is dependent on the environment at time \(t\)—as described by \(\epsilon(t)\)—as well as the vital rates that are realized by this environment and the current population density \(\boldsymbol{\theta}({\epsilon}(t), \mathbf{n}({t}))\). The matrix model determines the population at the next time step: 

\[{\textbf{n}}(t+1)={\textbf{A}}[{\boldsymbol{\theta}}(\epsilon(t), \mathbf{n}({t})),\epsilon(t)]{\textbf{n}}(t), \quad t=0,1,2,\ldots \tag1 \]

We use \(\epsilon\) to describe the proportional reductions in vital rates that result from an environmental disturbance.

We can define population recovery as occurring when the total population reaches a designated threshold \(N_{rec}\), such as the population’s carrying capacity or size prior to the disturbance. If we assume that a disturbance transpires at time \(t=0\), the recovery time will be the smallest integer solution to 

\[\prod_{i = 0}^{t} {\textbf{A}}[{\boldsymbol{\theta}}(\epsilon(i), \textbf{n}(i)),\epsilon(i)]{\textbf{n}}(0) \geq N_{rec}. \tag2 \]

The simplest way to describe a disturbance is with a step function, which dictates that the effects of a disturbance are either “on,” \(\epsilon(t) = \epsilon_0\), or “off,” \(\epsilon(t) = 0\). If we assume that the environment follows a step function and population growth is independent of density (as may be appropriate for endangered populations), the recovery time becomes the solution to the equation

\[\textbf{1}_m^\intercal\textbf{A}_0^{t-T_C}\textbf{A}^{T_C}_{\epsilon_0}\textbf{n}(0) = N_{rec}, \tag3 \]

where \(\textbf{1}_m\) is a \(m\times 1\) vector of ones.

How Sensitive are Recovery Predictions?

Sensitivity analysis measures the way in which small perturbations in a model parameter affect model output. When examining a population’s recovery, the recovery time’s sensitivity to a model input can help identify the most effective management or control strategy. One can calculate sensitivity in its simplest form just by taking a derivative. Using \((3)\) to model recovery, we can find the recovery time’s sensitivity by implicitly differentiating this equation. For instance, the sensitivity of the recovery time with respect to the magnitude of impact is given by 

\[\frac{dT_{rec}}{d \epsilon_0}=-\frac{\textbf{1}_m^\intercal\textbf{A}_0^{T_{rec}-T_C}(\textbf{n}(0)^\intercal\otimes \textbf{I}_m)\frac{d(\textbf{A}^{T_C}_{\epsilon_0})}{d\textbf{A}_{\epsilon_0}}\frac{d\mathrm{vec}[\textbf{A}_{\epsilon_0}]}{d\epsilon_0}}{(\textbf{A}^{T_C}_{\epsilon_0}\textbf{n}(0))^\intercal(\textbf{I}_m\otimes\textbf{1}_m^\intercal\textbf{A}_0^{-T_C})\frac{d\mathrm{vec}[\textbf{A}_0^{T_{rec}}]}{d{T_{rec}}}}, \tag4 \]

where the \(\otimes\) operator denotes the Kronecker product and the vec operator converts a matrix into a column vector by stacking the matrix’s columns. In a similar manner, we can also use \((3)\) to derive sensitivity formulas of the recovery time with respect to a vital rate or the initial population distribution [1].

Figure 1. Recovery time’s sensitivity to properties of the disturbance. 1a. Recovery time’s sensitivity to changes in survival reduction \(\epsilon_0\), assuming a 10-year duration of impact. 1b. Duration of impact \(T_C\), assuming a five percent reduction in survival. Figure adapted from [1].

To illustrate the utility of equations such as \((4)\), we present an application that investigates the recovery of sperm whales, which are impacted by a variety of disturbances that include oil spills and noise pollution. We use a discrete-time stage-structured model to examine a sperm whale population [3]. In Figure 1, we present the recovery time’s sensitivity for a sperm whale population with respect to changes in the magnitude of impact \(\epsilon_0\) and duration of impact \(T_C\). These graphs indicate that the recovery time is more sensitive to changes in \(\epsilon_0\) than \(T_C\). For instance, consider the effect of a 20 percent increase in these two parameters. If \(\epsilon_0=\) 0.20, a 20 percent increase in \(\epsilon_0\) heightens survival reductions from 20 to 24 percent, resulting in an additional 53 years of recovery time. In contrast, when we consider the recovery time’s sensitivity to changes in the duration of impact, this graph approaches a value close to six. Therefore, each additional year of impact increases the recovery time by approximately six years. If we thus increase \(T_C\) from 10 to 12 years (a 20 percent increase), the recovery time is only extended by 12 years.

Graphs such as Figure 1 provide important insights into a population’s recovery and management following a disturbance. In the context of sperm whales, Figure 1 suggests that conservation efforts should focus on reducing the magnitude of impact rather than the duration. In the case of a contaminant spill such as oil, this type of mitigation might include strategies that focus on removing the contaminant from the water. Similar analysis also indicates that management should concentrate on the mature individuals. As illustrated by this example, our model framework may provide useful insights even when limited data is available.

Azmy S. Ackleh is the dean of the College of Sciences and a professor in the Department of Mathematics at the University of Louisiana at Lafayette. He is also the Ray P. Authement Eminent Scholar and Endowed Chair in Computational Mathematics. Amy Veprauskas is an assistant professor in the Department of Mathematics at the University of Louisiana at Lafayette.