## Executive Summary

In the article by Keeling and Gilligan, Bubonic plague: A metapopulation model of a zoonosis, the authors begin by using a coupled system of SIR type models for humans, rats, and fleas to show how epidemics in the human population are preceded by an epidemic within the rat population. Keeling and Gilligan sought to have the most realistic model possible so they incorporated real data, and parameter sensitivity analysis when choosing the values for their parameters. They then performed an exhaustive analysis of this model, and conclude that it was adequate for short-term dynamics but something more robust would need to be developed for long-term analysis. Keeling and Gilligan then shifted their focus towards a stochastic model that was more biologically realistic. By combining the results of the model in question with modern biological data surrounding prevalence of common rodent species the risk of modern bubonic plague outbreaks were estimated probabilistically. All of this analysis provided them with a way to describe what a modern outbreak of the plague would look like (assuming that no measures were taken to combat the disease).

## Background

The bubonic plague is one of the most prolific and devastating diseases throughout global history. Spread by the bacterium Yersinia pestis, commonly found in flea populations, this disease has recurred for millenia. The first well documented epidemic was “Great Plague of Justinian”, occurring in the mid-sixth century this outbreak decimated Europe and Western Asia, spreading from Constantinople to what is now the United Kingdom. Likely the most well-known instance of the disease is the Black Death. Taking place in Europe during the 14th century this pandemic was responsible for the death of nearly half the population of Europe and lasted for over 100 years. A more contemporary outbreak occurred over large portions of Asia in the late 19th and early 20th centuries, this epidemic persisted in India for about 50 years, claiming over 10 million lives before eventually arriving in San Francisco via stowaways on shipping routes.

These devastating outbreaks, and the knowledge that disease persists at a low level today show the need to have accurate and available models to determine not just the spread of the disease but, more importantly, the risks to humans and the ways the disease can outbreak into human populations. With the disease models readily available the persistence of Yersinia pestis can be modeled, but to gain a high degree of accuracy in relation to historically seen attributes, multiple models must be extended and combined.

The SIR model originated with the Kermack-McKendrick model put forth in the early twentieth century as continuous compartment model with susceptible, infective, and removed classes. This model provides the framework for the expanded system derived by Keeling and Gilligan, where they use the typical classes for both rat and human populations, as well creating two classes for flea populations: average number of fleas living on an individual rat, and number of fleas searching for a host – a major driving force for the disease. As aforementioned, the termination of rats drives the fleas to find new hosts, often ending up on humans or pets, thereby transmitting the disease to a new susceptible population.

## The Model

Keeling and Gilligan developed an SIR like model to describe the disease dynamics within the human population, within the rat population and within the flea population. The results of their development are shown below:

**Brief Description of the terms in the model**

System for the human population:

First equation - Susceptible Humans:

- Susceptible humans are born at a rate proportional to the number of susceptible and recovering humans. Susceptible humans die at a rate proportional to the number of proportional humans.
- Susceptible humans become infected at a rate proportional to the number of susceptible humans times the number of free searching fleas times the "Force of the infection" (detailed below).

Second Equation - Infective Humans:

- The first term is identical to the last term in the previous equation.
- The infected humans leave either because they die at a rate proportional to the number of infected humans or they move to the recovered class after a period of time at a rate proportional to the number of infected humans.

Third Equation - Recovered Humans:

- A fraction of the humans leave the infected class and enter the recovered class after having been sick for a period of time at a rate proportional to the number of infected humans. It should be noted that this fraction who enter the recovered class is the survival rate.

System for the Rat population:

First Equation - Susceptible Rats:

- Rats are born at a rate proportional to the number of susceptible rats, but this growth rate saturates as the rat population approaches the carrying capacity for the rats.
- Susceptible rats are also born at a rate proportional to the number of recovering rats times the probability that the infant rats did not inherit immunity to the disease from their recovering parents.
- The susceptible rates die at a rate proportional to the number of susceptible rats.
- The susceptible rates leave the susceptible class and join the infected class at a rate proportional to the number of susceptible rats divided by the total number of rats times the force of the infection.

Second Equation - Infective Rats:

- The first term is the same as the last term in the previous equation.
- The rats leave the infected class either because they die at a rate proportional to the number of rats or they leave the infected class after having been sick for the average time a rat is sick.

Third Equation - Recovered Rats:

- The removed rats are born at a rate proportional to the number of removed rats who give birth to offspring that are immune to the disease.
- The infected rats move to the removed class at a rate proportional to the number of rats that have been infected for the average infectious period.
- The removed rats die at a rate proportional to the number of removed rats.

System for the flea population:

First Equation - Average Number of Fleas on a Rat:

- The flea index (average number of fleas living on a rat) increases at a rate proportional to the flea index, but saturates as the value approaches the carrying capacity for the fleas.
- The last term represents free fleas finding a host, and adding to the flea index.

Second Equation - Infective Fleas Searching for a Host:

- The number of free fleas increases at a rate proportional to the number of infected rats times the flea index. It should be noted that this proportion is the natural death rate of the rats times the rate that the infected rats die. The number of free fleas decrease as the fleas die at a rate proportional to the number of free fleas.

**Force of Infection**

Keeling and Gilligan noted that the human population dynamics do not have an affect on the other two populations so they instead focused on the number of fleas who failed to find a host in conjunction with the total rat population. They argued that these free, searching fleas could feed on and subsequently infect a human host, or pass the disease to other animals cohabiting with humans. They then derive the value for the potential force of infection (it is apparent that the true force of infection would be less because not all fleas will actually infect a human). This term is given below:

## Determining Parameter Values

According to Keeling and Gilligan most of the parameters were taken from
their literature, experiments and observations.
The parameters that they did not have information on were chosen to be within
what they called, “biologically reasonable bounds”.
These values are listed in the following table:

**Sensitivity of Parameters**

Keeling and Gilligan then carried out a sensitivity analysis on these parameters. Namely they tracked the change in three model outputs as the parameter values were varied. These three outputs were the equilibrium number of rats, the period of the first epizootic cycle, and the potential force for that first cycle. This measure S was defined as follows:

They then used this relationship in calculating the sensitivities of the parameters, and their results are summarized in the following chart:

## Deterministic Results

The motivation for developing this model was to determine how the rat and flea population could lead to an outbreak of the plague in the human population, and the result is summarized in the following plot:

We can see that the force of infection for a given number of free fleas changes exponentially as we vary the number of rats, and that for a given number of rats the force of infection changes linearly as we change the number of free fleas. We can also observe that the force of infection is greatest when there are a lot of free searching fleas, but only a handful of rat hosts. This result is logical, as a large searching flea population with few rats to find as hosts have a greater potential to find human hosts.

The force of infection and total number of infective rats were plotted as a function of time by Keeling and Gilligan. This plot is analogous to findings from the computational analysis of the purported model, where the force of infection was found and plotted as a function of time. See the model below to the left, and the found results below to the right (note :both plots are logarithmic in the y-axis).

Keeling and Gilligan noticed that the force of infection lagged behind the number of rat cases by four weeks, and this agreed with the literature that they had been referencing. They also noticed that the human epidemics were damped much faster than the rodent epidemics, and this means that for there to be human cases there must have been a large epidemic in the rodent population.

The deterministic model can accurately describe the short term behavior of the plague, but as can be seen in the plot we only observe fixed point behavior. This is an obvious flaw because throughout history there have been seemingly spontaneous outbreaks of the plague, and this model would never predict such outbreaks. This motivated Keeling and Gilligan to seek another description that would allow eruptions of the disease after the disease becomes endemic in the rat population.

**Student Evaluation of the model**

To evaluate this model a set of initial conditions needed to be determined. In order to gain relevance and realism a scaled version of populations from Boulder County, Colorado were used. Taking population data from the 2014 U.S. Census Bureau's population estimates and using disease-carrying rodent densities of North America from a study by Gary Witmer and Gilbert Proulx approximations for local human and rodent populations were formulated and scaled. Given a rodent density of 300 per hectare and the population estimates for Boulder the initial population sizes become 313,000 humans, and 1,880,000 rodents. Based on flea to rat ratios taken from Keeling and Gilligan a reasonable estimate of 5,000,000 free fleas was used. These numbers were then reduced by a factor of 1000, giving 313 humans, 1880 rats, and 5000 fleas. A scenario was chosen in which only one of the scaled "humans" would carry the disease, and 20% of the rat population was infective initially. These initial conditions can be seen summarized in the following table:

A fourth order Runge-Kutta method was then used to numerically evaluate the system of differential equations, allowing plots to be formulated for the human and rat populations, as well as for the force of infection over time, which is shown in the section above. The plots for the human and rat populations over time are shown below, where the populations are broken into susceptible, infective, and recovered classes.

Here it is worth noting that there is an initial epidemic within the rats that causes an epidemic in the human population. Then there is a second epidemic in the rats that does not cause an epidemic in the humans, this is likely caused by the smaller rat population during the second epizootic wave . This is exactly what Keeling and Gilligan predicted would happen, and it supports the idea that there is a need for a model in which multiple epidemics can happen within the human population. It should also be noted that the size of subsequent epidemics in the rat population will continue to decrease, and this further supports the observation by Keeling and Gilligan that this model exhibits fixed point behavior.

## Greater Biological Realism

Standard assumptions have historically been that plague outbreaks were triggered by a specific event, be that infected rats transferred via ships, or an infected human entering a local population, which could accurately be shown by the deterministic model described above. The following model puts forth an alternative scenario, one which has been seen in past outbreaks, where the bubonic plague can erupt randomly. Keeling and Gilligan hypothesize that outbreaks such as these arise from movements in the rat population, where the disease is endemic.

To gain increased realism the authors introduce three forms of heterogeneity: temporal forcing, individuality, and spatially distinct subpopulations. The latter of these allows the total population to be discretized and coupled via random movements of rats and fleas. The stochasticity of the model now allows for random fluctuations in the system further explaining much of the historically seen behaviors of the bubonic plague.

One important result from this extension is that outbreaks of the plague do not necessarily need to arise from a specific triggering event as previously assumed. In larger populations the disease can become endemic in the rat population, posing little to no threat to humans, but following the right series of random events, the force of infection could rise and another wave of human cases could erupt. Namely, following a surge in the susceptible rat population the disease has the capability of jumping between species and posing serious health risks to the human population, which has been seen historically, but is not effectively shown by the deterministic model.

## Analysis/Interpretation

The overall purpose of Keeling and Gilligan’s work was to provide an accurate model in order to not only qualify previous episodes of bubonic plague outbreak but to evaluate risk for potential outbreaks in current human populations. Through sensitivity analysis and the introduction of stochasticity into the model one parameter combination was determined to be of supreme importance when assessing the risk to human health from bubonic plague, this is the term aKr where a is the searching efficiency of fleas, and Kr is the rat (or more generally: rodent) carrying capacity. Since the searching efficiency and the carrying capacity should move inversely to one another as the size of a subpopulation changes their product, aKr will remain (close to) constant. What this leads to, interpretively, is a scenario where a large population of susceptible rats in a subpopulation can become infective rapidly, where the mass expiration of rats in the population leaves infective fleas to search for new hosts thereby transmitting the disease to new populations.

Once spread to human populations the disease moves through the population via the differential equations shown above with a force determined by the total number of rats in the system, and the number of free infected fleas. As shown in the stochastic model this can arise from random fluctuations in the endemic state of bubonic plague among the rat population, and does not need to be driven by an outside force. This is a significant improvement to the deterministic model, as the introduction of stochasticity shows behaviors seen historically but not shown by the original model, namely the spontaneous eruption of disease in isolated populations of susceptible humans.

## References

[1] Keeling, M. J., and C. A. Gilligan. "Bubonic Plague: A Metapopulation Model of a Zoonosis." Proceedings of the Royal Society B: Biological Sciences 267.1458 (2000): 2219-230. Web.

[2] Witmer, Gary and Gilbert Proulx. Rodent outbreaks in North America. (2010. Pages 253-267 in Grant R. Singleton, Steve R. Belmain, Peter R. Brown, and Bill Hardy, editors. Rodent outbreaks: ecology and impacts. Los Banos (Philippines): International Ride Research Institute.

[3] Annual Estimates of the Resident Population: April 1, 2010 to July 1, 2014 Source: U.S. Census Bureau, Population Division Release Dates: For the United States, regions, divisions, states, and Puerto Rico Commonwealth, December 2014. For counties, municipios, metropolitan statistical areas, micropolitan statistical areas, metropolitan divisions, and combined statistical areas, March 2015. For Cities and Towns (Incorporated Places and Minor Civil Divisions), May 2015.

[4]Ligon, B. Lee. "Plague: A Review of Its History and Potential as a Biological Weapon." Seminars in Pediatric Infectious Diseases 17.3 (2006): 161-70. Web.