Project 12.5.3 - Flu Epidemic

For this problem in our final project we will be looking at the spread of a disease through a population that is fixed in size with no persons entering or leaving the population. To solve this problem we will be using a continuous SIR model. In this model we split the population into three different groups, the susceptible population (S), the infected population (I), and the removed population (R). We start with every member of the population being either susceptible of being infected or infected. As members of the population are infected they can then spread the flu to others and once a person is removed (already had the flu or death) they can not get the infection again.

To build the model for our problem we need the following information:

We can now build our coupled SIR model for the spread of our disease. \[ \frac{dS}{dt} = -0.002 I(t) S(t) \] \[ \frac{dI}{dt} = -0.7 I(t) + 0.002 I(t) S(t) \] \[ \frac{dR}{dt} = 0.7 I(t) \]

Euler’s Method

To solve this model we will use Euler’s method for solving differential equations with a step size of 0.5. In Euler’s method we end up up using the information from our initial conditions and our rate of change formulas to estimate the values for susceptible, infected, and removed for the next time step. The basic formulation for any Euler’s method problem is. \[ x_n = x_{n-1} + f(t_{n-1}, x_{n-1}, y_{n-1}) * \Delta t \] \[ y_n = y_{n-1} + g(t_{n-1}, x_{n-1}, y_{n-1}) * \Delta t \] Where \(f(t, x, y)\) and \(g(t, x, y)\) are the rate of change functions for \(x\) and \(y\) and \(\Delta t\) is our step size.

In our coupled SIR model we get the following equations that define our values at the current time step:
\[ r_n = r_{n-1} + (0.7 * i_{n-1}) \Delta t \] \[ i_n = i_{n-1} + (-0.7 * i_{n-1} + 0.002 * i_{n-1} * s_{n-1}) \Delta t \] \[ s_n = s_{n-1} + (-0.002 * i_{n-1} * s_{n-1}) \Delta t \]
With \(\Delta t = 0.5\) or half of a week.

We have chosen to implement this solution in the following R function, flu(), that steps through the solution. We have also created two helper functions. The first, graph(), creates a ggplot of the model run. The second, tbl(), cleans up the model and outputs the model run as a table for printing.

library(ggplot2)
library(reshape2)
library(pander)
library(dplyr)
flu <- function(R, I, S, a = 0.7, b = 0.002, delta_t = 0.5, itt = 24){
  # Computes the output of the SIR Model for Epidemics. The total population = S+I+R
  #
  # Args:
  #   R: Removed, the number of members of the population that have are no longer susceptiable to the disease.
  #   I: Infected, the number of members of the population that have the disease.
  #   S: Susceptiable, the number of members of the population that can become infected.
  #   a: Removal rate per week, defaults to a value of 0.7
  #   b: Transmision coefficient, defaults to a value of 0.002
  #   delta_t: The step size, defaults to 0.5
  #   itt: The number of itterations in weeks to run the model, defaults to 24 weeks.
  #
  # Return:
  #   Outputs the table of outputs from the model run and a plot of the three populations over the time period.
  week <- seq(0, itt, by = delta_t)
  removed <- c()
  susceptiable <- c()
  infected <- c()
  for(i in 0:(itt/delta_t)){
    R_tmp <- R
    I_tmp <- I
    S_tmp <- S
    removed <- c(removed, R)
    susceptiable <- c(susceptiable, S)
    infected <- c(infected, I)
    R <- R_tmp + (a*I_tmp)*delta_t
    I <- I_tmp + (-a*I_tmp + b*I_tmp*S_tmp)*delta_t
    S <- S_tmp + (-b*I_tmp*S_tmp)*delta_t
  }
  df <- data.frame(week, removed, infected, susceptiable)
  return(df)
}

graph <- function(df){
  # Helper function for the flu function to graph the output of the model run.
  #
  # Args:
  #   df: the data frame output from the flu model run.
  #
  # Return:
  #   ggplot of the model run.
  melted <- melt(df, id.vars = "week")
  ggplot(melted, aes(x = week, y = value, color = variable)) + geom_point() + geom_line()
}

tbl <- function(df){
  # Helper function for the flu function to output a cleaned table of the model run.
  #
  # Args:
  #   df: the data frame output from the flu model run.
  #
  # Return:
  #   pander table of the output from the flu model run, the output is cleaned to only return 
  #     whole weeks.
  clean <- df %>%
    filter(week %% 1 == 0)
  panderOptions("digits", 4)
  pander(clean)
}

First Run - I(0) = 2

For our first run of the model we will be looking at a run with initial conditions of \(I(0) = 2\) and \(S(0) = 998\).

flu_df <- flu(R=0, I = 2, S=998, a = .7, b = 0.002)
tbl(flu_df)
week removed infected susceptiable
0 0 2 998
1 1.854 5.425 992.7
2 6.872 14.59 978.5
3 20.3 38.36 941.3
4 55.09 94.94 850
5 138.2 202.1 659.7
6 301.6 311.4 387
7 523.6 296 180.4
8 713.2 191 95.81
9 829.9 103.6 66.47
10 892.2 52.69 55.16
11 923.6 26.09 50.31
12 939.1 12.77 48.1
13 946.7 6.219 47.07
14 950.4 3.02 46.57
15 952.2 1.465 46.33
16 953.1 0.7104 46.22
17 953.5 0.3443 46.16
18 953.7 0.1669 46.13
19 953.8 0.08086 46.12
20 953.8 0.03919 46.12
21 953.9 0.01899 46.11
22 953.9 0.009201 46.11
23 953.9 0.004459 46.11
24 953.9 0.00216 46.11 
graph(flu_df)

We see that when we start with 2 people infected with the flu we get a peak number of infections around week 6 with approximately 323 people infected by the flu. This run of the model ends with 954 people infected and 46 people who remain susceptible to the illness. We reach a level of less than one person infected by week 16 and the flu has run its course by that time.

Second Run - I(0) = 5

For our second run of the model we will be looking at a run with initial conditions of \(I(0) = 5\) and \(S(0) = 995\).

flu_df <- flu(R=0, I = 5, S=995, a = .7, b = 0.002)
tbl(flu_df)
week removed infected susceptiable
0 0 5 995
1 4.629 13.49 981.9
2 17.05 35.63 947.3
3 49.44 88.99 861.6
4 127.7 193 679.3
5 285 307.4 407.6
6 506.4 303.1 190.5
7 701.7 199.4 98.94
8 823.7 108.9 67.38
9 889.2 55.47 55.35
10 922.3 27.48 50.23
11 938.6 13.45 47.91
12 946.6 6.544 46.83
13 950.5 3.176 46.31
14 952.4 1.54 46.06
15 953.3 0.7458 45.94
16 953.8 0.3612 45.88
17 954 0.1749 45.85
18 954.1 0.08469 45.84
19 954.1 0.04101 45.83
20 954.2 0.01985 45.83
21 954.2 0.009613 45.83
22 954.2 0.004654 45.83
23 954.2 0.002254 45.83
24 954.2 0.001091 45.83 
graph(flu_df)

We see that when we start with 5 people infected with the flu we end up getting slightly different results then when we started with just two people infected. We reach a maximum number of people infected at one time of 325 between during week five. This slightly increases the total number of people that are infected at one time and also moves the peak of the number of people a week earlier then on our first run. We end up with 955 people infected and 45 that are still susceptible. The flu also runs it course by week 15 when we drop to less then one person infected, this is also a week earlier then we saw in the previous run of the model.

Third Run - I(0) = 10

For our third run of the model we will be looking at a run with initial conditions of \(I(0) = 10\) and \(S(0) = 990\).

flu_df <- flu(R=0, I = 10, S=990, a = .7, b = 0.002)
tbl(flu_df)
week removed infected susceptiable
0 0 10 990
1 9.24 26.73 964
2 33.7 68.53 897.8
3 94.81 157.6 747.5
4 227.1 281.9 491
5 438.4 322.5 239.1
6 651.6 232.8 115.6
7 795.5 131.7 72.85
8 874.9 67.88 57.24
9 915.4 33.77 50.79
10 935.5 16.54 47.91
11 945.4 8.049 46.58
12 950.2 3.903 45.94
13 952.5 1.89 45.64
14 953.6 0.9145 45.49
15 954.1 0.4423 45.42
16 954.4 0.2139 45.39
17 954.5 0.1034 45.37
18 954.6 0.05001 45.36
19 954.6 0.02418 45.36
20 954.6 0.01169 45.36
21 954.6 0.005654 45.36
22 954.6 0.002734 45.36
23 954.6 0.001322 45.36
24 954.6 0.0006391 45.36 
graph(flu_df)

We again see some shifting of the model when we increase the number of people infected at the start of the run to 10. The max number of people infected is still around 325 but now this happens near the end of week 4. We still end up with 955 people infected and 45 that are still susceptible. The flu now runs it’s course with less than one person infected by week 14 which is also a reduction of one week from the previous model run.

Fourth Run - I(0) = 100

For our last run of the model we will be looking at a run with initial conditions of \(I(0) = 100\) and \(S(0) = 900\)

flu_df <- flu(R=0, I = 100, S=900, a = .7, b = 0.002)
tbl(flu_df)
week removed infected susceptiable
0 0 100 900
1 89.25 226.3 684.5
2 274.2 356.2 369.6
3 525.9 322.5 151.5
4 729.3 194.6 76.12
5 846.9 100.5 52.65
6 906.7 49.24 44.02
7 935.9 23.64 40.42
8 949.9 11.26 38.82
9 956.6 5.337 38.08
10 959.7 2.526 37.74
11 961.2 1.195 37.58
12 961.9 0.5647 37.5
13 962.3 0.2669 37.47
14 962.4 0.1261 37.45
15 962.5 0.05962 37.44
16 962.5 0.02817 37.44
17 962.5 0.01331 37.44
18 962.6 0.006292 37.44
19 962.6 0.002973 37.44
20 962.6 0.001405 37.44
21 962.6 0.000664 37.44
22 962.6 0.0003138 37.44
23 962.6 0.0001483 37.44
24 962.6 7.008e-05 37.44 
graph(flu_df)

Increasing our initial condition to 100 people infected has caused some interesting changes in the model. Our first and most noticeable is that we have increased the total number of people infected at once to 363 people and this now happens in the second week of the infection. We also end up with 963 people infected and 37 people still susceptible to the disease. The flu has also run its course by week 12 which is a two week reduction from the previous run of the model.

Improved Euler’s Method

A second methodology for solving this type of problem is to use an use The Improved Euler’s method. In this method we use the same ideas of estimating the slope and using that to predict the next value, however we average the current function value with the new function value to improve the estimate.
\[ x_{n+1}^* = x_{n} + f(t_{n}, x_{n}, y_{n}) * \Delta t \] \[ y_{n+1}^* = y_{n} + g(t_{n}, x_{n}, y_{n}) * \Delta t \] \[ t_{n+1} = t_n + \Delta t \] \[ x_{n+1} = x_n + [f(t_{n}, x_{n}, y_{n}) + f(t_{n+1}, x_{n+1}^*, y_{n+1}^*)] * \frac{\Delta t}{2} \] \[ y_{n+1} = y_n + [g(t_{n}, x_{n}, y_{n}) + g(t_{n+1}, x_{n+1}^*, y_{n+1}^*)] * \frac{\Delta t}{2} \]

We will modify our original flu() function to use the Improved Euler’s method for solving this type of problem.

imp_flu <- function(R, I, S, a = 0.7, b = 0.002, delta_t = 0.5, itt = 24){
  # Computes the output of the SIR Model for Epidemics. The total population = S+I+R
  # This version uses the Improved Euler's Method for solving.
  # Args:
  #   R: Removed, the number of members of the population that have are no longer
  #       susceptiable to the disease.
  #   I: Infected, the number of members of the population that have the disease.
  #   S: Susceptiable, the number of members of the population that can become infected.
  #   a: Removal rate per week, defaults to a value of 0.7
  #   b: Transmision coefficient, defaults to a value of 0.002
  #   delta_t: The step size, defaults to 0.5
  #   itt: The number of itterations in weeks to run the model, defaults to 24 weeks.
  #
  # Return:
  #   Outputs the table of outputs from the model run and a plot of the three populations 
  #       over the time period.
  week <- seq(0, itt, by = delta_t)
  removed <- c()
  susceptiable <- c()
  infected <- c()
  for(i in 0:(itt/delta_t)){
    R_tmp <- R
    I_tmp <- I
    S_tmp <- S
    removed <- c(removed, R)
    susceptiable <- c(susceptiable, S)
    infected <- c(infected, I)
    R_star <- R_tmp + (a*I_tmp)*delta_t
    I_star <- I_tmp + (-a*I_tmp + b*I_tmp*S_tmp)*delta_t
    S_star <- S_tmp + (-b*I_tmp*S_tmp)*delta_t
    R <- R_tmp + ((a*I_tmp) + (a*I_star))*(delta_t/2)
    I <- I_tmp + ((-a*I_tmp + b*I_tmp*S_tmp) + (-a*I_star + b*I_star*S_star))*(delta_t/2)
    S <- S_tmp + ((-b*I_tmp*S_tmp) + (-b*I_star*S_star))*(delta_t/2)
  }
  df <- data.frame(week, removed, infected, susceptiable)
  return(df)
}

First Run - I(0) = 2

For our first run of the model we will be looking at a run with initial conditions of \(I(0) = 2\) and \(S(0) = 998\).

impflu_df <- imp_flu(R=0, I = 2, S=998, a = .7, b = 0.002)
tbl(impflu_df)
week removed infected susceptiable
0 0 2 998
1 2.646 6.87 990.5
2 11.64 23.02 965.3
3 40.79 71.08 888.1
4 122.4 174.1 703.5
5 284.5 271 444.6
6 481.1 263.1 255.8
7 643.8 194.8 161.5
8 755.2 127.2 117.6
9 825.4 78.39 96.24
10 867.8 46.92 85.24
11 893 27.66 79.32
12 907.8 16.17 76.05
13 916.4 9.407 74.2
14 921.4 5.459 73.15
15 924.3 3.163 72.54
16 926 1.831 72.2
17 926.9 1.06 72
18 927.5 0.6131 71.88
19 927.8 0.3546 71.81
20 928 0.2051 71.77
21 928.1 0.1186 71.75
22 928.2 0.06859 71.74
23 928.2 0.03966 71.73
24 928.2 0.02294 71.73 
graph(impflu_df)

We see that when we start with 2 people infected with the flu in our improved model we get results that are slightly different to what got with our original model. We get a peak number of infections around week 5 with approximately 280 people infected by the flu. This run of the model ends with 928 people infected and 72 people who remain susceptible to the illness. We reach a level of less than one person infected by week 18. We note that this is a lower peak and a more spread out infectious period then what we got using the regular Euler’s method.

Second Run - I(0) = 5

For our second run of the model we will be looking at a run with initial conditions of \(I(0) = 5\) and \(S(0) = 995\).

impflu_df <- imp_flu(R=0, I = 5, S=995, a = .7, b = 0.002)
tbl(impflu_df)
week removed infected susceptiable
0 0 5 995
1 6.583 16.96 976.5
2 28.37 54.14 917.5
3 92.81 144.1 763.1
4 235.7 255.7 508.5
5 432.2 275.1 292.6
6 607.6 214 178.4
7 731.7 142.9 125.4
8 811.1 88.96 99.97
9 859.4 53.51 87.07
10 888.2 31.62 80.2
11 905.1 18.5 76.42
12 914.9 10.77 74.3
13 920.7 6.251 73.1
14 924 3.622 72.41
15 925.9 2.096 72.01
16 927 1.212 71.78
17 927.7 0.7011 71.65
18 928 0.4053 71.57
19 928.2 0.2343 71.53
20 928.4 0.1354 71.5
21 928.4 0.07829 71.49
22 928.5 0.04525 71.48
23 928.5 0.02616 71.48
24 928.5 0.01512 71.47 
graph(impflu_df)

Once again we see that increasing the number of people infected by the flu at the start of the model run moves the peak of the illness earlier and reduces the total run of the disease. We still maintain a peak number of infections of 280 during week 4 and a finish with 929 people infected and 71 people still susceptible to the flu. We reach less than one person infected in week 17.

Third Run - I(0) = 10

For our third run of the model we will be looking at a run with initial conditions of \(I(0) = 10\) and \(S(0) = 990\).

impflu_df <- imp_flu(R=0, I = 10, S=990, a = .7, b = 0.002)
tbl(impflu_df)
week removed infected susceptiable
0 0 10 990
1 13.06 33.23 953.7
2 54.44 98.28 847.3
3 161.7 214.3 624
4 345.4 282.9 371.7
5 538 245.8 216.2
6 685.1 172.4 142.5
7 782.2 109.8 108.1
8 842.2 66.74 91.03
9 878.2 39.64 82.13
10 899.4 23.26 77.3
11 911.8 13.55 74.61
12 919 7.867 73.09
13 923.2 4.557 72.22
14 925.6 2.636 71.72
15 927 1.524 71.44
16 927.8 0.8806 71.27
17 928.3 0.5087 71.18
18 928.6 0.2939 71.12
19 928.7 0.1697 71.09
20 928.8 0.09803 71.07
21 928.9 0.05662 71.06
22 928.9 0.0327 71.06
23 928.9 0.01889 71.05
24 928.9 0.01091 71.05 
graph(impflu_df)

We again see some shifting of the model when we increase the number of people infected at the start of the run to 10. The max number of people infected is still around 283 but now this happens near the start of week 4. We still end up with 929 people infected and 71 that are still susceptible. The flu now runs it’s course with less than one person infected by week 16 which is also a reduction of one week from the previous model run.

Fourth Run - I(0) = 100

For our last run of the model we will be looking at a run with initial conditions of \(I(0) = 100\) and \(S(0) = 900\)

impflu_df <- imp_flu(R=0, I = 100, S=900, a = .7, b = 0.002)
tbl(impflu_df)
week removed infected susceptiable
0 0 100 900
1 114.4 236 649.6
2 319.7 315.5 364.7
3 532.3 266.7 201.1
4 689.4 181.7 128.8
5 790.5 112.9 96.62
6 851.6 67.29 81.14
7 887.5 39.26 73.22
8 908.3 22.66 68.99
9 920.3 13 66.67
10 927.2 7.438 65.37
11 931.1 4.246 64.64
12 933.3 2.422 64.23
13 934.6 1.38 64
14 935.3 0.7865 63.87
15 935.8 0.448 63.79
16 936 0.2552 63.75
17 936.1 0.1454 63.72
18 936.2 0.0828 63.71
19 936.3 0.04716 63.7
20 936.3 0.02686 63.7
21 936.3 0.0153 63.69
22 936.3 0.008712 63.69
23 936.3 0.004962 63.69
24 936.3 0.002826 63.69 
graph(impflu_df)

Increasing our initial condition to 100 people infected has caused some interesting changes in the improved Euler’s model. Our first and most noticeable is that we have increased the total number of people infected at once to 315 people and this now happens in the second week of the infection. We also end up with 936 people infected and 64 people still susceptible to the disease. The flu has also run its course by week 14 which is a two week reduction from the previous run of the model.

Conclusions

We see from the output of both of the models that the spread of the flu in the fixed population follows a consistent pattern. We start off with a rapid increase in the number of people infected over the first few weeks. We reach a peak number of infections between weeks 4 and 6, and the flu runs it’s course in 12 to 18 weeks. The only runs which showed deviations from this pattern were runs where the number of people infected started at 100 individuals. In both models this caused the infection to reach a higher peak number of infected and caused more of the population to be infected by the end of model run.

Comparing the two models we see that we get very similar results when using Euler’s Method and the Improved Euler’s Method. Since the Improved Euler’s method is designed to give more accurate results we would assume that in this situation our first model gave us more of a worst case scenario with more people infected at once and more people infected by the disease at the end of the model run. In a health care scenario we could easily argue that this is actually the better version of the model to run. Often in this setting we need to build systems and infrastructure to handle the worst case plus some and our Euler’s Method model gives us that. The advantage to the results of the Improved Euler’s method model, is that it spreads the infection out over a longer period. We could argue that we may want to use this model to determine the length of time that the flu is active . By combining the outputs of both of these models we get a picture for the worst case of what could happen with the flu infection in a given year. This was a very interesting model to put together and we learned some very interesting things about the way that coupled models work in a mathematical modeling.