SIR Model of the Black Plague in 14th Century Europe

Co-Authors: Katie Ericson, Julia Goyco, Renae Rhode, Lauren Smith

Advisor: Emek Kose

Mathematical Modeling - Fall 2019

Abstract

The devastation of the Black Plague in 14th century Europe is well-known. We know how many people died and how long the plague lasted. However, what if humans never learned how to save themselves? What if humans never figured out how to stop the plague from spreading? How long would it have taken for the Black Plague to completely destroy Europe? This paper answers these questions by utilizing general SIR models to successfully describe the Black Plague’s effect on all of Europe, if people had never learned the basic principle of cleanliness. After establishing the assumptions, like birth rate being the same as natural death rate and using only human to human interactions, the equations were run through a differential equation solving program to map the spread of the disease to all people in Europe. The results of this modelling showed a dramatic change in length and severity from the historical plague to our modelled plague.

Introduction

Everyone that has ever taken a World History class has heard about the Black Death. This epidemic spread throughout all of Europe in the mid 1300’s and impacted history in an unforgettable way. The Black Death is known to stem from the plague, specifically the bubonic plague, and is caused by an infection from the bacterium Yersinia pestis (Howard, 2019). The plague is thought to have been brought from China to Europe and spread dramatically via lice, rats, and humans.

The symptoms of the Black Death include vomiting, nausea, fever, and growths on the skin called buboes ("Black Death | Causes, Facts, and Consequences", 2019). The key to surviving this plague was as easy as washing hands, washing clothes, and popping the buboes. However, in the 1300’s, people didn’t understand germs. They thought that praying the illness away would be the cure, but unfortunately it wasn’t. From 1347 to 1351 ("Black Death | Causes, Facts, and Consequences", 2019), the Black Plague devastated Europe, leaving almost nothing but infected corpses behind to continue to infect others. 

As the Black Plague ran through all of Europe, people slowly began learning what not to do. The two most important discoveries were popping the buboes and burning the clothes of the deceased. As previously stated, people didn’t understand germs. Indicating that when they wore the clothing of the deceased, they were becoming infected as well. Once, these people understood that concept, infectivity rate decreased. It was widely known that confidence in religion dropped dramatically as religion wasn’t saving everyone from the plague and the economy of Europe was decreasing in prosperity over the duration of the disease ("Black Death | Causes, Facts, and Consequences", 2019).

The main goal of this paper is to explore the devastation the Black Death caused Europe during the 1300’s and to see the impact of never learning how to survive the cure. This paper will model the spread of the Black Death through Europe and its impact on the population, assuming that the populace never learned to clean their clothes or wash their hands. Seeing how quickly the plague could have wiped out all of Europe creates an interesting comparison with actual history and how long the plague lasted in real life. 

Previous Works

Most previous works about modeling the Black Plague describe the interaction between the plague and the rodents/fleas involved in the process. The paper “Human ectoparasites and the spread of plague in Europe during the Second Pandemic” (Dean, Krauer & Walloe, 2018) uses a compartmental model to illustrate how the plague is transmitted through fleas and body lice, a compartmental model for human-to-human transmission and a compartmental model for rats transmitting the plague amongst themselves and other humans (Dean, et al., 2018). We incorporated several of the parameters - such as the human-to-human contact rate, the average low infection period, and average high infection period - from this previous work in our model. 

Although we focused on modeling the plague in Europe in the 1300s, there are other papers that study different outbreaks, “Epidemiological analysis of the Eyam plague outbreak of 1665-1666,” which used an SEIR model that accounted for a human-to-human infection rate for people living in the same household as an infected.

In a 2000 paper “Bubonic plague: a metapopulation model of a zoonosis,” the authors treat the modeling of the black plague as a zoonosis model rather than a compartmental model. A zoonosis model describes a disease that is passed from animal to humans (like rabies and ebola) which was used for the plague because transmission from rats to humans (via fleas) created the largest amount of infected (Zoonosis).

The black plague is an important disease to model because it can be applied to future models of diseases as well as help us understand a historically important disease that killed millions over many years in Europe as well as outbreaks in other continents.

The Model

Model 1

Figure 1: Compartmental Model 1

We started off with a basic SIR model. A couple of changes we made were that we had two paths leading from susceptible to infected, as well as two paths leading from infected to removed. This is due to people getting infected by both living (ꞵ ) and dead (ƛ) individuals who have (and had) the plague. These, however, have different rates, because some people were more likely to come into contact with the dead (having to bury them). The two paths leading from infected to removed represent the two different infection periods. The majority of the literature averaged the infection rates, but we found that there were generally two different types of infection: short time frame (𝛔) and long time frame (𝛄). The short time frame was about 2 days, and the long time frame was about 8 days, so we made them two different rates.

Model 2

Figure 2: Compartmental Model 2

In order to more accurately predict the untempered spread of the plague, we introduced two other variables, ⍺ and κ. The first variable that we introduced, ⍺, was our capping variable. ⍺ was the initial population of Europe at the time, 85 million, divided by two. This didn’t come from any single source, but was more of a generalization for us to use in our modeling to try and reach a more accurate result. To fit into our R code better, we created a second new variable, κ, that was just ⍺ (not changing) plus our S variable (changes with model). In our equations for the model, we divided both ꞵ  and ƛ by κ. This is due to ꞵ  and κ being our two different infection rates. 

Assumptions

Results

Model 1

SIR Model 1

Figure 3: Model 1

Initially, our model did not cap the infection range a single infected individual had. This would basically be saying that each individual in the susceptible category would be equally likely to come into contact with an infected individual. This meant that the number of susceptibles dropped extremely quickly, due to the initial number of susceptibles being so large. Similar to the zombie model that we discussed in class, if you start off with a small population, this is not as much of an issue because it mirrors the real-life interactions more closely. However, if you are considering all of Europe as the population, this does not really work out. Since the infections happen so quickly, having more susceptibles just means it will be spread faster. The initial model can be seen in Figure 3, with the t-axis being time in weeks, and the y-axis being number of individuals.

Model 2

SIR Model 2

Figure 4: Model 2

This next model takes a completely different shape than the previous model. In fact, the susceptible curve is almost completely inverted from what our first model showed. In the first model, the susceptibles dropped quickly, then leveled out at zero; whereas, the revised model shows the susceptibles remaining generally constant until week 55, or so, and then going over a “roller coaster” type hill that lasts about 20 weeks. If our graph was extended on the x-axis, the susceptibles leveling out to zero in the current model would also be shown. These changes are due to the additional parameters that were included in the system of differential equations in our attempt to increase the accuracy of the result. 

As shown in Figure 4, capping the infection rate made a much more accurate model. Peak infection occurred at about week 90 from the time it was introduced. You can see that once a certain number were infected, there is a dramatic decrease of susceptibles due to the epidemic. By week 100 (about two years), all of Europe had been killed due to the plague. In each model, we end up with a “doomsday” equilibrium, where everyone dies. The main difference is that it just takes different amounts of time, but from two weeks to roughly two years is a major difference.

Model 3

Figure 5: Model 3

Modeling doesn’t always go as planned and there are many opportunities for learning experiences. What this group has learned is that it is very important for multiple people to run the modeling code and check each other’s work. This is the model that has been affectionately named the “Why?” model, shown in the figure below. Oddly enough, the susceptible curve appears to be the same shape in each model. The equations and parameters should be the same for both of the models, but apparently, something was amiss that has yet to be found. What this model is essentially telling us is that a spike in infections near week one cause a spike in deaths that do not affect the number of susceptibles until week 4700, or so. By week 5000, all of Europe has died. This is approximately 96 years. We still are not sure if R was short-circuiting, or if we input any of the information wrong, but what we do know is that the model is wrong. 

Sensativity Analysis for Model 1 and Model 2

Sensitivity analysis is used to decide the strength of model outcomes to parameter values since there are errors in assumed parameter values. For the sensitivity analysis, we fixed the values of the parameters and observed how the entire system is affected by incrementally changing only one of the parameters. This analysis provides an estimate for how much the parameters will affect the system and are therefore their sensitivity. 

To test the sensitivity of each parameter to the different variables we started with a fixed susceptible population with one infected person and zero recovered people. All of the parameters were fixed at the values in the table while varying only the one being tested by plus and minus ten percent. Then the S, I, and R equations were graphed. The graphs were analyzed to see if the parameter had any effect on the susceptible, infected and recovered populations. 

We found that ƛ and ⍺ are sensitive parameters, they affect the outcome of the model the most. When ƛ is increased by 10% it takes about 90 weeks for the entire population to end up in the recovered box. The susceptible population decreased at a faster rate while the infected population increased faster and there was a slight increase in the number of infected individuals. When ƛ is decreased 10% it takes more than 100 weeks for the entire population to end up in the recovered box. When ƛ is smaller it takes, more time form the disease to spread. This makes sense because the more bodies there are to come into contact with, the more it will spread. When ⍺ was increased, the infected population increases at a slower rate. The disease takes longer to spread when ⍺ is increased. When ⍺ is decreased the disease spreads more quickly. This is an expected outcome because  is part of the capping expression (⍺ + S), the smaller it gets, the more people interact and spread the disease.

Perturbing 𝛔 and 𝛄 gives a slight change in the outcome of the model. This means that these parameters are not very sensitive. There was no visible change in the model output when ꞵ was increased or decreased by 10%.  This was an unexpected result since ꞵ is the human to human infection rate so if it increased the disease should spread more quickly.

Stability Analysis

The equilibrium points were found by setting the equations to zero and using Maple software to solve for the points. The stability of the equilibrium points was analyzed by using the Jacobian method. In this method, the eigenvalues and eigenvectors of the Jacobian matrix of the system were computed using Maple. The equilibrium points for both models are (0, 0, 0), (0, 0, 85,000,000), (85,000,000, 0, 0). The stability of the equilibrium point is determined by looking at the eigenvalues. If the eigenvalues are all positive (sink) or all negative (source) then the point is stable. If the eigenvalues are all negative (source) or if they are positive and negative values (saddle) then the point is unstable. The eigenvectors show the direction in which the system is increasing or decreasing. 

Stability for Model 1

The Jacobian matrix for Model 1 is given by:


Substituting  S = 0, I = 0 and R = 0 the eigenvalues of the matrix are 0, 0 and -1.39. The eigenvectors are (1, 0, 0), (0, 0, 1) and (0, -1, 1). The equilibrium point is unstable and is decreasing in the direction (0, -1, 1).

Substituting S = 0, I = 0 and R = 85,000,000 the eigenvalues are -29,750,000, -1.402, and 0.0000000005608. The eigenvectors are (0,0,1), (0, -1, 1) and (2.14 x 10⁷, -2.14 x 10⁷,1). This equilibrium point is unstable and is decreasing in the (0,0,1)  and (0,-1,1) direction and increasing in the (2.14 x 10⁷, -2.14 x 10⁷,1) direction.

Substituting S = 85,000,000, I = 0 and R = 0 the eigenvalues are 0, 4.24 x 10⁶, and -9,73. The point is unstable.

Stability for Model 2

The Jacobian matrix for Model 2 is given by:


Substituting  S = 0, I = 0 and R = 0 the eigenvalues of the matrix are 0,0 and -1.39. The eigenvectors are (1, 0, 0), (0, 0, 1) and (0, -1, 1). The equilibrium point is unstable and is decreasing in the direction (0,-1,1).

Substituting S = 0, I = 0 and R = 85,000,000 the eigenvalues are -1.39, -0.233 and 5.102⨉10⁻¹¹. The eigenvectors are (0, 0, 1), (-0.832, -0.168, 1) and (0, -1, 1). This equilibrium point is unstable and is decreasing in the (0,0,1)  and (-0.832, -0.168,1) direction and increasing in the (0, -1, 1) direction.

Substituting S = 85,000,000, I = 0 and R = 0 the eigenvalues are 0, -1.56, and 0.207. The eigenvectors are (1, 0, 0), (-1.51, 0.15, 1) and (0.13, -1.13, 1) The point is unstable and is decreasing in the (-1.51, 0.15,1) direction and increasing in (0.13, -1.13, 1) direction.

Conclusion

In this paper, we examined the spread of the black plague throughout 14th century Europe. Using a simple SIR model to map the epidemic, we accounted for the passage of the disease through human-to-human and body-to-human contact in order to maintain a straight-forward and easy to understand approach to the question. 

Our initial model did not produce a realistic result, in that the entirety of the population of Europe was infected and died from the disease within a couple of hours. In order to account for the fact that not all individuals that were susceptible would be equally likely to come into contact with an infected individual, we introduced a capping variable in our second model to produce a much more realistic spread. 

From our sensitivity analysis, we discovered that ƛ (body-to-human) and ⍺ (capping variable) are the most sensitive. When ƛ increases it takes less time for the disease to spread, which is an understandable result: the more bodies there are to come into contact with, the more it will spread. When ⍺ is decreased, the disease again would spread at a faster rate than without it. This also makes sense as ⍺, a part of the capping expression (⍺ + S), is a fraction, the smaller it gets, the more people interact and spread the disease. From the stability analysis, both Model 1 and Model 2 were found to be unstable at all equilibrium points. 

Our model demonstrates the doomsday equilibrium, which occurs at about 90 weeks or 2 years, and not the more accurate representation of a period of prolonged illness of the population for 4 years before recovery. This is due to the fact that we assumed there would be no entry to the population, there is no immunity (and therefore, people who could not get sick), and that people would not learn to adapt and figure out the “cure” and prevention to the disease. We also did not take into account separate infection rates for the lice and rat population, assuming that the spread of the disease would be constant and not vary with the corresponding population of infectors. This meant that the disease spread faster and was more “effective” as it infected additional members of the susceptible population. 

In conclusion, our SIR model is an effective method of modeling a simplified version of the Black Plague as it occurred in 14th century Europe. The model also demonstrates that without the variation in human biology as well as the ability to adapt and find solutions to various ailments, the population of Europe would have quickly been wiped out.

Future Works

The models introduced in this paper are the simplified expressions of a complicated disease. Of the restrictions for the models, limiting the interactions (of infecteds) and means of infection were the most necessary in order to prevent an increase in complexity. We examined only the idea of infection passing through contact with the bodies of the deceased and human-to-human. Since we wanted to exclude the other means of spreading and obtaining the disease, we ignored the lice and rat populations of Europe and their interactions with the infected and susceptible populations. Given more time, we could research the rat (living and dead) and lice populations of 14th century Europe to account for the spread of the disease through those methods.

One of the reasons the disease eventually died out in Europe (without killing the entire population) is the discovery of the impact of popping the cysts (buboes) and washing hands in order to cure and impede the spread of the disease. In order to include these in the model, we would need to find the approximate time at which these discoveries were made in order to map their effect on the disease, as they would not be immediate and would take effect gradually and after a period of unimpeded infection. 

Another factor that could be taken into consideration for future models would be the movement of Europeans during that time period (1347-1351), as it we assumed a mostly static population. We would also factor in natural birth and death rates (as we presumed they were equal in order to focus on the impact of the plague deaths.

References