Section 20.1 Models of an epidemic
A basic model of epidemic divides the population into three groups: susceptible \(S\) (able to be infected), infected \(I\text{,}\) and removed \(R\) (not able to be infected). Often group \(R\) is called “recovered” but it also contains deceased and vaccinated. One can modify the model in various ways: include a separate group for deceased, allow for a possibility of reinfection, etc. The sum \(S+I+R\) remains constant in time, because other demographic factors such as births do not substantially change the picture of an epidemic.
It is convenient to express \(S, I, R\) as a proportion of total population. This means these quantities are between \(0\) and \(1\) and their sum is \(1\text{.}\) An important assumption of the model is that the rate of infection is proportional to the product \(SI\text{.}\) This is based on the simplifying idea that any two people in the population are equally likely to meet and transmit the disease, and since \(SI\) is the number of all susceptible-infected pairings, the number of transmissions will be proportional to it. This is not too far from reality when “population” means, for example, student population on some campus. When modeling the population of a large country, one would have to consider the effects of geography.
Thus, the rate of new infections is \(\beta SI\) where \(\beta\) is the coefficient of proportionality that depends on how contagious the disease is and what preventive measures are taken by the population. At the same time, the rate of new recoveries can be reasonably assumed to be \(\gamma I\text{,}\) proportional to the number of currently infected. So we have the system
with some initial conditions \((S_0, I_0, R_0)\) where typically \(R_0=0\) and \(I_0\) is small. The constants \(\beta,\gamma\) are the parameters of this system. Their values have a major impact on the long-term behavior of solutions; yet these values are very difficult to estimate in practice (especially \(\beta\)). One might try to estimate them by studying a small-scale outbreak closely and fitting a model to the observed counts of infected and recovered people.
We could solve the system (20.1.1)-(20.1.2)-(20.1.3) numerically using the methods of Chapter 19 but to keep the focus on the modeling aspect, we use Matlab's solver ode45. Its basic syntax is
[t, y] = ode45(rhs, [a, b], y0);
where rhs is the right hand side of the ODE system, [a, b] is the interval on which solution must be found, and y0 is the initial condition at the left endpoint of that interval. The initial condition can be a single number (when solving a single equation) or a vector (when solving a system). It is preferable to use a row vector for y0 when solving a system, because of how the solution is presented by ode45: see below. Note that rhs must take two variables in this order: rhs = @(t, y) ... even if the system does not involve time (is autonomous). The output of function rhs must be a column vector (it can be a single number if we are solving one equation).
The solver returns a vector of \(t\)-values and a vector or matrix of corresponding \(y\)-values. It is a matrix if we solve a system; in this case y(:,1) approximates the first unknown function, y(:,2) the second, and so on. The \(t\)-values are not uniformly spaced in general; they are chosen based on the behavior of the system. One can plot the solution with plot(t, y).
Sometimes it is desirable to compute the solution at specific points of interval \([a, b]\text{,}\) for example at equally spaced points with step \(h\text{.}\) This is achieved by replacing the interval parameter with a vector of \(t\)-values, for example
h = 0.1; [t, y] = ode45(rhs, a:h:b, y0);
Example 20.1.1. Implement SIR model.
Use ode45 to solve the equations (20.1.1)-(20.1.2)-(20.1.3) on the interval \([0, 100]\) with \(\beta = 0.3\text{,}\) \(\gamma = 0.2\text{,}\) and \(\mathbf y_0 = (0.99, 0.01, 0)\text{.}\)
beta = 0.3;
gamma = 0.2;
y0 = [0.99, 0.01, 0];
rhs = @(t, y) [-beta*y(1)*y(2); beta*y(1)*y(2) - gamma*y(2); gamma*y(2)];
[t, y] = ode45(rhs, [0, 100], y0);
plot(t, y, 'LineWidth', 3)
legend('Susceptible', 'Infected', 'Recovered')
In this scenario a part of the population (less then half) escapes the infection because the epidemic ends on its own after the transmission rate drops. Try varying the coefficients \(\beta, \gamma\) and observe the effect.
In reality, \(\beta\) is not necessarily constant in time. Let us introduce seasonality into the model by making \(\beta\) a periodic function of \(t\text{,}\) for example
beta = @(t) 0.3 + 0.2*cos(0.2*t);
which is a periodic function that varies between \(0.1\) and \(0.5\text{.}\) Try this out (note that beta should be replaced by beta(t) in the formula for rhs).