Math Insight

1. If a person with influenza recovers with a probability of $\gamma = 0.2 \text{ day}^{-1}$, on average, how many days will a person be infective? Expression your answer in terms of $\gamma$ and also give the number of days.

   There are a total of $N=100,000$ people, which mean that we could have, at most, $100,000$ susceptible people. In order to give an upper bound on how fast the infection will spread, let's imagine for the moment that we have $100,000$ susceptibles. (Even though the number of susceptibles must be less than $100,000$ as soon as we have one infective, the number of susceptibles could still be close to $100,000$ when there are just a few infectives.) For each of these $100,000$ susceptible people, there is a probability $\beta = 4.0 \cdot 10^{-6} \text{ day}^{-1}$ that they will contract influenza from a single infective. During the course of those days when that person is infective, how many people, on average, will contract influenza from that infective person? This number is called the basic reproductive number, and we denote it $R_0$. Give a number for $R_0$ based on the given values of the parameters. Also, give a general formula for $R_0$ in terms of $N$, $\beta$, and $\gamma$.

   Is $R_0$ larger or smaller than 1? Based on the value of $R_0$, do you expect that the disease will spread or die out?

2. The total number of individuals $N = S + I + R$ should be constant. Verify that, for the model you derived, it is indeed true that $\diff{N}{t} = \diff{S}{t} + \diff{I}{t} + \diff{R}{t} = 0$.
3. Since $N$ is constant, we really don't need to keep track of all three variables $S$, $I$, and $R$. If we know what $S$ and $I$ are, what must $R$ be? (Solve the equation defining $N$ for $R$.) Since we don't need to keep track of $R$ separately, we can ignore the $\diff{R}{t}$ equation and reduce the system to just the two equations for $\diff{S}{t}$ and $\diff{I}{t}$. Write the new system of two equations. The new equations shouldn't depend on $R$, just on $S$ and $I$ (and the total number $N$). Replace any occurrences of $R$ with the expression you determined for $R$. (Again, write the equations in terms of the parameters $\beta$, $\gamma$, and $\alpha$.)
4. Calculate the $I$-nullcline, that is the equation where you set the rate of change of $I$, $\diff{I}{t}$, to zero. You can substitute the values of $N$, $\beta$, and $\gamma$ if you like, so that everything is numbers except for $S$ and $I$.

   Factor this equation, so you have the product of two factors equaling 0. Since at least one of those factors must be zero, you have two possible ways for $\diff{I}{t}$ to be zero. One condition should involve only $S$ and the other only $I$, so they will be vertical and horizontal lines on the $S-I$ phase plane. Plot both pieces of the $I$ nullcline on the phase plane.
5. Calculate the $S$-nullcline, i.e., the equation where $\diff{S}{t}=0$. Solve this equation for $I$. The result will be a function of $S$. (It will be a fraction with an expression involving $S$ in both the numerator and denominator.) This function will not be simple to sketch by hand on the above phase plane. Use a computer program or a graphing calculator to determine the shape of the function, and sketch the $S$-nullcline (use a different color or line style or label it to distinguish it from the $I$-nullcline) on the above phase plane. Pay careful attention to the locations where the $S$-nullcline intersects the $I$-nullcline. It should intersect in two places.
   How to plot a function in R (Show)(Hide)

   How to plot a function in R

   You don't need to use R to plot the $S$-nullcline. But, if you want to use R, here's how you can plot the function. First, define $S$ to be the vector of all the values of $S$ you want to plot. Looking at the above axes, it seems as though we want to go up to $S=110,000$. In R, to be make $S$ be a vector of points from 0 to 110,000, counting by 100s, you can this command.
   S = seq(0, 110000, 100)

   Next, you want to make $I$ be a vector with the corresponding values. For example, if the $S$-nullcline were $I=\frac{0.5S^3}{1+S^2}$ (which it isn't), you could create the vector $I$ with the following command.
   I=0.5*S^3/(1+S^2)

   With these vectors in place, you use use the plot command, using the vector S for the first argument (the values of $x$) and the vector $I$ for the second argument (the values of $y$). The simplest command would be
   plot(S,I,type='l')
   but that would make the $I$ axis be a much larger scale than on the above phase plane. To limit $I$ between 0 and 20000, you can use the ylim, as $I$ is the $y$-axis.
   plot(S,I,type='l',ylim=c(0,20000))

   (Hide)

6. Next, calculate the equilibria, i.e., the points $(S_e,I_e)$ where both $\diff{S}{t}$ and $\diff{I}{t}=0$. These are the points where the $S$-nullcline intersects the $I$-nullcline.

   From the phase plane graph, you can see that the $S$-nullcline intersects each piece of the $I$-nullcline once. To determine the coordinates of each of these equilibria, plug in the value of $S$ or $I$ from the $I$-nullcline equation into the $S$-nullcline equation. Then, solve for the other variable. Draw and label the equilibrium points on the above phase plane.

7. The nullclines divide the biologically reasonable portion (with $S \ge 0$ and $I \ge 0$) of the phase plane into four regions. For each region, pick any point $(S,I)$ in that region. Calculate the rates of change, $\diff{S}{t}$ and $\diff{I}{t}$, at each of those points. (Really all you need to determine is the sign of each rate of change.) Draw those points on the phase plane and sketch a direction vector from those points. Each direction vector should start at the point $(S,I)$ and indicate the direction $(\diff{S}{t}, \diff{I}{t})$, which is the direction the solution should move when crossing that point. The magnitude of the vector is irrelevant; just make sure that it points in the approximate correct direction (such as north-east, north, south-west, etc.).
8. The $I$-nullcline cuts the (biologically realistic portion of) the $S$-nullcline into two pieces. Pick a point $(S,I)$ on each piece and calculate direction vectors for each. Plot the points and direction vectors on the phase plane.
9. The $S$-nullcline cuts the vertical portion of the $I$-nullcline into two pieces. Pick a point $(S,I)$ on each piece and calculate direction vectors for each. Plot the points and direction vectors on the phase plane.

   The $S$-nullcline cuts the horizontal portion of the $I$-nullcline into two pieces, but only one piece has realistic values of $S$ that are less than or equal to $N$. For that piece, pick a point, calculate a direction vector, and plot the results on the phase plane.

10. Calculate the Jacobian matrix of the system. Evaluate the Jacobian matrix at each equilibrium and calculate the eigenvalues. Use the result to classify each equilibrium.
11. Use this information to sketch the solution trajectory $(S(t),I(t))$ on the phase plane for the initial condition with 100 people infective and the rest susceptible, i.e., for $(S(0),I(0)) = (99900, 100)$. In each of the four regions of the phase plane, the solution should move in the approximate direction of the direction vector you sketched and it should cross the nullcline segments parallel to their corresponding direction vector. If a solution gets near an equilibrium, the behavior of the solution should then match the classification of that equilibrium.

12. Given the solution you sketched in the phase plane, sketch plausible plots of $S(t)$ and $I(t)$ versus time. You don't have enough information to estimate the timing, so the scale on the $t$-axis is arbitrary. However, make sure the relationship between $S(t)$ and $I(t)$ qualitatively match your solution in the phase plane. Pay particular attention to when the variables are increasing or decreasing. Their long term behavior should match the long term behavior in the phase plane.

Your goal is to eradicate influenza from this population of 100,000 individuals. You have a number of intervention strategies available to you.

1. You could give yearly vaccinations to 60% of the population, keeping them immune to influenza at all times. For the model, this intervention would multiply the population size $N$ by 0.4 to change it to $N=40,000$.
2. You could reduce the rate of infection by 60% by educating the population in good hygiene (such as washing hands and covering coughs and sneezes) and encouraging sick people to avoid contact with others. For the model, this intervention would multiply the infection rate parameter $\beta$ by 0.4 to change it to $\beta=1.6 \cdot 10^{-6}$.
3. You could administer drugs to help the infectives fight off the infection, increasing the rate of recovery by 150%. For the model, this intervention would multiply the recovery rate parameter $\gamma$ by 2.5 (equivalently, divide it by 0.4) to change it to $\gamma=0.5$.
4. You could administer drugs that decrease the rate of loss of immunity to influenza by 60%. For the model, this intervention would multiply the immunity loss parameter $\alpha$ by 0.4 to change it to $\alpha=0.004$.

Using the following steps, analyze the effectiveness of each strategy.

1. For each of these strategies, calculate the basic reproduction number $R_0$. Compare these values to the original $R_0$ and predict the success of the interventions from this value.

2. In the original model, you should have determined that there is a disease-free equilibrium $(S_e,I_e) = (N,0)$, which corresponds to everyone being healthy, albeit susceptible. However, that equilibrium was unstable so that the introduction of even a small number of infectives lead to a flu outbreak. Once that outbreak started, the state of the system progressed toward a different equilibrium, one where the disease persisted in the population.

   Pick one intervention that your $R_0$ calculations predict would be successful and one intervention that these calculations predict would be unsuccessful. For each of these interventions, classify the disease-free equilibrium. From this classification, determine whether or not the introduction of a small number of infectives should lead to an outbreak of the flu.

   Also, for both interventions, determine if there is another biologically plausible equilibrium (meaning $S_e \ge 0$, $I_e \ge 0$, and $S_e + I_e \le N$). If so, determine the stability of that equilibrium.

   Based on these results, predict whether or not the intervention should succeed in eradicating the flu.

3. For your chosen intervention that should eradicate the flu, sketch the phase plane, along with the nullclines, equilibria, and sample direction vectors. Sketch some solutions on the phase plane, using initial conditions with a few representative values of $I(0)$ and setting $S(0)=N-I(0)$ (i.e., start with no one in the removed class). For all initial conditions, do these solutions tend to the disease free equilibrium as time increases? Evaluate the success of the intervention.