An SIS model with periodic contact rate

Let's explore what happens when the contact rate $\beta$ is a periodic function of time. For example, for a childhood disease we might imagine that $\beta$ is higher during the school year and lower during the summer. Following the paper ``Asymptotic behavior in a deterministic model'' by H.W. Hethcote, Bull. Math. Biol. 35:607-614, 1973, let's take

$$\beta(t) = 2 - 1.8 \cos (5 t).$$ (4)

We thus want to solve

$$\frac{dI}{dt} = (\beta(t) N - \alpha) I - \beta(t) I^2.$$ (5)

This is known as a nonautonomous differential equation because the righthand side has explicit time dependence. It can be rewritten as an autonomous system of equations (in which the righthand side has no explicit time dependence) by defining a new independent variable $\tau$ satisfying $d\tau/dt = 1$. Our system of equations is then

$$\frac{dI}{dt} = (\beta(\tau) N - \alpha) I - \beta(\tau) I^2,$$ (6)

$$\frac{d\tau}{dt} = 1.$$ (7)

This system of equations may be solved numerically with the following Matlab programs. First SIS.m:

global N alpha

alpha = 1;
N = 1;

options = odeset('MaxStep',0.01);
[T,Y] = ode45('func_SIS',[0 10],[0.2 0],options);  

figure(1)
plot(T,Y(:,1));
xlabel('t');
ylabel('I');

for i=1:1000         %also plot a scaled version of beta
   tt(i) = 0.01*i;
   bb(i) = beta(tt(i));
end
hold on;
plot(tt,0.1*bb,'r');

Text version of this program

Next, func_SIS.m:

function dy = func_SIS(t,y)

global N alpha

I = y(1);
tau = y(2);

dI = (beta(tau)*N - alpha)*I - beta(tau)*I^2;
dtau = 1;

dy = [dI;dtau];

Text version of this program

Finally, beta.m:

function r = beta(t)

     r = 2 - 1.8*cos(5*t);

Text version of this program

Figure 1 shows the numerical solution for this system with $N=1$ and $\alpha =4$ in blue, and $0.1*\beta (t)$ in red. It is seen that $I$ decays to zero with a few wiggles in the transient. Figure 2 shows the numerical solution with $N=1$ and $\alpha =1$, again with $0.1*\beta (t)$ in red. Here $I$ settles into a periodically oscillating state with period equal to the period of $\beta$; note that the maximum of $I$ does not occur at the same time as the maximum of $\beta$.

Figure 1: Numerical solution with $N=1$, $\alpha =4$ in blue, $0.1*\beta (t)$ in red.

Figure 2: Numerical solution with $N=1$, $\alpha =1$ in blue, $0.1*\beta (t)$ in red.

To understand these results, Hethcote defines the average reproduction number as

$$\bar{R}_0 = \frac{\bar{\beta} N}{\alpha},$$ (8)

where the bar indicates the time average. If $\bar{R}_0 \le 1$ then $I$ damps in an oscillatory manner to zero, whereas if $\bar{R}_0 > 0$ then $I$ approaches a periodic solution. For Figure 1, $\bar{R}_0 = 0.5$, and for Figure 2, $\bar{R}_0 = 2$, so we find the expected behavior.