An SEIR model with periodic contact rate

Following the paper by Aron and Schwartz, let's assume that

$$\beta(t) = \beta_0 (1 + \beta_1 \cos 2 \pi t).$$ (14)

Furthermore, let's take

$$\mu = 0.02, \qquad \alpha = 35.84,$$

$$\gamma = 100, \qquad \beta_0 = 1800.$$

The following Matlab programs generate numerical solutions for this system. First SEIR.m:

global alpha beta0 beta1 gamma mu

mu = 0.02;
alpha = 35.84;
gamma = 100;
beta0 = 1800;
beta1 = 0.2;

options = odeset('MaxStep',0.01);
[T,Y] = ode45('func_SEIR',[0 50],[0.0658 0.0007 0.0002 0],options);  

figure(1)
plot(T,-log(Y(:,3)));
xlabel('t');
ylabel('-ln(I)');

figure(2);
plot(-log(Y(:,1)),-log(Y(:,3)));
xlabel('-ln(S)');
ylabel('-ln(I)');

Text version of this program

Next, func_SEIR.m:

function dy = func_SIS(t,y)

global alpha beta0 beta1 gamma mu

S = y(1);
E = y(2);
I = y(3);
tt = y(4);

dS = mu - beta_SEIR(tt)*S*I - mu*S;
dE = beta_SEIR(tt)*S*I - (mu + alpha)*E;
dI = alpha*E - (mu+gamma)*I;
ds = 1;

dy = [dS;dE;dI;ds];

Text version of this program

Finally, beta_SEIR.m:

function r = beta_SEIR(t)

   global alpha beta0 beta1 gamma mu

   r = beta0*(1+beta1*cos(2*pi*t));

Text version of this program

Playing around with the initial conditions and different values of $\beta_1$, one finds the following solutions. Note that these figures show the solution after integrating for a long time - therefore, they give the asymptotically stable behavior. Also, it is most convenient to plot minus the logarithm times the variables.

Figure 3: Numerical solution in phase space for $\beta _1=0.05$.

Figure 4: Time series for the numerical solution for $\beta _1=0.05$.

Figure 5: Numerical solution in phase space for $\beta _1=0.2$.

Figure 6: Time series for the numerical solution for $\beta _1=0.2$.

Figure 7: Numerical solution in phase space for $\beta _1=0.26$.

Figure 8: Time series for the numerical solution for $\beta _1=0.26$.

Figure 9: Numerical solution in phase space for $\beta _1=0.27$.

Figure 10: Time series for the numerical solution for $\beta _1=0.27$.

We see that at $\beta _1=0.05$, the period of $I$ is the same as the period of $\beta$, namely 1. At $\beta _1=0.2$, the period of $I$ is twice the period of $\beta$. A period doubling bifurcation has occurred between these values of $\beta_1$. Increasing $\beta_1$ to 0.26, another period doubling bifurcation occurs, and at $\beta _1=0.26$ the period of $I$ is four times the period of $\beta$. Between $\beta _1=0.26$ and $\beta _1=0.27$, an infinite number of such period doubling bifurcations occur (this is called a period doubling cascade). At $\beta _1=0.27$ the solution is chaotic: it never exactly repeats itself!

A very important property of chaotic systems is that there is sensitive dependence on initial conditions, i.e., solutions which are close to each other at a given time will eventually diverge away from each other. Try integrating this model for two nearby initial conditions to see this property. In fact, you're encouraged to spend some time playing around with this model - you'll find a lot of surprises!