Exemplary Final Term Paper 6

Author: Renato de Oliveira

Introduction

Modelling Epidemic Dynamics

With the advent and employment of vaccination procedures, many diseases that posed as great factors of mortality (especially infant mortality) started being controlled. In order to improve the capacity of combating the morbidity caused by major epidemics that affect human populations, quantitative studies of population dynamics have been extensively developed, gathering as much of data and theory as possible.

For this purpose, quantitative descriptions of epidemic dynamics have brought highly relevant understanding of the key features of disease outbreak, spread and disappearance. Not only on the field of epidemics, but many other under the life sciences, mathematical modelling has provided deep insights. On epidemic dynamics particularly, theories on disease modelling were central to controlling diseases such as polio and smallpox .

Measles, in particular, has a historic that provides for rich and reliable epidemiological databases, thus making it a paradigm on the study of epidemics. It is a well characterized disease, with a great impact on society and on public health through the time, even as of today, as it still causes mortality in developing countries. These features make measles a well suitable choice for a standard for modelling epidemics . Moreover, the failure of mass vaccination campaigns to eradicate it in developig countries allows for the continued interest in researching this disease and its impact on society.

On the article under study , the authors explore the properties of an established mathematical model of epidemic dynamics, the SIR model. It comprises a compartmentalization of the population (and takes its name from it) in discrete fractions, namely the Susceptible, whose individuals may be subject to contagion, the Infective, whose individuals carry the disease and are able to transmit the epidemics, and the Recovered, whose individuals are healed from the infective population and don't become infective again. The SIR model provides easily approachable tools to understand disease dynamics while still accounting satisfactorily for accurate predictions on a basic level. There exists a variety of altered or updated versions of this model to satisfy more realistically demanding or conditionally constrained assays, whereas the basic SIR model itself allows for very relevant and impactant conclusions upon its essential framework (this will be discussed on further sections).

The Authors' Approach

The authors make an extensive consideration about the mathematical aspects of the SIR model both untouched and under adaptations to vaccination policies. They chose to evaluate the model theoretically for the Constant Vaccination strategy, which is a well established immunization protocol worldwide, and for the Pulse Vaccination strategy

Pulse Vaccination Strategy

The Pulse Vaccination strategy comprises periodic vaccinations in the population, not only on newborns, but on all age cohorts. The same ratio of the population is vaccinated between established periods of time. This impulsive vaccination is thought to interfere directly with the dynamics of the epidemic as it oscillates near extinction. This given, Pulse Vaccination policy might prove advantageous, more cost-efficient than Constant Vaccination policy.

Hypotheses

In order to evaluate the above propositions, the authors guide their efforts by the following hypotheses:

• Through pulsed vaccination policy, a smaller proportion of the population needs to be immunized; This would make the pulse vaccination strategy more efficient than constant vaccination alone.
• If the time between vaccination pulses is smaller than some threshold, the epidemic will reach an acceptable stable state or disappear; The dynamics of the SIR model under no combative strategies or under constant vaccination strategy predict that the infectious population stabilizes above zero, thus permitting posterior outbreaks of disease. It would then be great to prove that the pulse vaccination strategy can eventually lead to eradication of disease.

Reproducing Results

In our paper we build working implementations of the system that are flexible enough to permit alterations and adaptations, thanks to the power of Mathematica. We implement the system under the three main conditions proposed in the article, unaffected SIR model, SIR under constant vaccination strategy and SIR under pulse vaccination strategy, as well as implement additions to the dynamics of the system, such as seasonal enforcing and mixing both strategies of vaccination.

We also aid graphically in the description of the system under these various conditions, using the standard parameters for measles they provide in the article. Each major proposition is confirmed by simulations we ran and plotted as our results.

Generally speaking we were able to demonstrate the results shown in the article, expand its description visually and extend over its analysis of disease dynamics. These will be detailed in the remaining sections.

Model Description

The Pure SIR Model

The SIR model assesses the evolution of a population under an epidemic with time. It compartmentalizes the population into three groups, each assigned a state variable. It thus assumes that these three strict groups describe every state in which a individual from the population can be found at.

As the groups represent fractions of the population, all three state variables sum up to one, so we can discard one of them arguing that we can extract its information from the difference of the total and the other state variables. This approach effectively reduces the model's number of dimensions to two. For this purpose, we thus assume that a recovered individual will not become infected again.

The state variables are $S(t)$, the proportion of susceptible population, $I(t)$, the proportion of infectious population and $R(t)$, the proportion of population that recovered from the infectious state (and will not become infected again).

The equations of state contain the three parameters described in Table 1. The state variables have the following rate equations, detailed in table 2.

(1.1)	$\frac{dS}{dt} = m - (\beta{}I + m)S$
(1.2)	$\frac{dI}{dt} = \beta{}IS - (m+g)I$
(1.3)	$\frac{dR}{dt} = gI-mR$

The system is effectively reduced to two dimensions by observing the relation $R(t)+S(t) + I(t) = 1$. Another key assumption of the model is that the total birth and total death rates ($m$) are the same.

Table 1. Main parameters used in state equations.

Symbol	Name	Standard value for measles	Description
$m$	Rate of birth	0.02	The rate at which individuals are born (which equals the rate of death).
$\beta$	Contact rate	1800	The rate at which infectious individuals actually contaminate susceptibles.
$g$	Rate of recovery	100	Inverse of the mean infectious period, the rate at which infectious individuals recover from illness (and acquire immunity).

Table 2. Detailing the differential equations governing the model.

Equation	Term	Meaning
(1.1)	$-(\beta{}I + m)S$	$\beta{}I$ determines the capacity of the infectious population to contaminate susceptibles, and this capacity is increased as the susceptible population itself increases, while the $m$ term represents the rate at which the susceptible population dies, thus both multiply $S$.
	$m$	Represents the rate at which susceptible individuals are born, thus increasing the size of the susceptible population.
(1.2)	$\beta{}IS$	This term was present in equation (1.1) subtracting from $m$, and represents the rise (or decrease) in the infectious population.
	$-(m+g)I$	The death rate of the entire population ($m$) and the rate of recovery of infected individuals ($g$) both decrease the infectious population.
(1.3)	$gI$	This is the rate of recovery of infected individuals by the total number of infectious individuals, representing the direct growth of recovered population.
	$-mR$	This is the total death rate by the total recovered population, representing the direct decrease in the number of recovered individuals.

The pure, unaffected SIR model has two equilibrium points, one trivial, infetion-free equilibrium point, and a non-trivial, epidemic equilibrium point. We set the rate equations (1) to zero and solved for the state variables. This yields the coordinates of the system at which it remains in a stable state, as equaling all the rates of change to zero means the state variables are not going to be altered, thus keeping the system at the same state. Letting $S$ be represented along the abcissae and $I$ be represented along the ordinates, we found the following coordinates for the equilibrium points

(2) $(1,0),(S^{*},I^{*})^{ }_{ }$, where $S^{*}$ and $I^{*}$ are treated as parameters in order to facilitate code notation. They are detailed in Table 3.

Table 3. Composed parameters used in coding.

Symbol	Name	Formula	Standard value for measles	Description
$R_{0}$	Basic reproductive rate of disease	$\frac{1}{S^{*}} \text{ } \text{or} \text{ } \frac{\beta}{m+g}$	17.9964	A conjugation of the three basic parameters $m$, $g$ and $\beta$, which corresponds to the inverse of $S^{*}$, the "epidemic" equilibrium value of $S$.
$S^{*}$	Epidemic equilibrium value of $S(t)$.	$\frac{m + g}{\beta}$	0.0555667	The value the state variable $S(t)$ for the equilibrium point of the epidemic state, i.e. when the proportion of infectious population $I(t)$ is greater than zero.
$I^{*}$	Epidemic equilibrium value of $I(t)$.	$\frac{m(R_{0} - 1)}{\beta}$	0.000188849	The value the state variable $I(t)$ assumes at the equilibrium point of the epidemic state, i.e. when its value is greater than zero.

Reproducing Results

In order to reproduce the analysis and explore the main features of the pure SIR model, we generate phase plots, nullclines and trajectories of the system making use of the standard measles parameters.

We use Mathematica's annotation power to solve (through Solve) for the equilibrium points and set them as parameters to facilite further coding. In order to produce nullclines, we use the functions Solve and Reduce to solve the rate equations for each state variable when their rates of change are set to zero. Using VectorPlot we generate a slope field and plot the nullclines, equilibrium points and a trajectory (generated by ParametricPlot) to graphically exploit the properties of the system.

The Constant Vaccination Strategy

The constant vaccination strategy consists of continually vaccinating a proportion $p$ of the population as they are born. This is reflected in the equations of state as a decrease on the rate of birth of the susceptible population, $S(t)$, as $pm$ newborns will immediately integrate the recovered population. The rate of change of the susceptible population then becomes

(3) $\frac{dS}{dt} = (1-p) m - (\beta I + m) S$.

The stability of the equilibrium points and their values are shifted according to the parameter p. There is a critical value of p (equation (4)) for which stability changes according to equations (5) and (6).

(4)	$p_{c} = 1 - \frac{1}{R_{0}}$
(5)	$S^{*}_{0} = (1-p)$, $I^{*}_{0}=0$,	for $p > p_{c}$.	$(S^{*}_{0},I^{*}_{0})$ is stable.
(6)	$S^{*'}_{1}=S^{*}_{1}$, $I^{*'}_{1}=I^{*}_{1} - \frac{m}{m+g} p$,	for $p < p_{c}$.	$(S^{*'}_{1},I^{*'}_{1})$ is stable.

$(S^{*}_{0},I^{*}_{0})$ is the infection-free equilibrium and $(S^{*'}_{1},I^{*'}_{1})$ is the epidemic equilibrium.

Reproducing Results

Once again making use of Mathematica's great annotation power we use Solve and Reduce to search for the new equilibrium points. We find them to satisfy equations (5) and (6).

Next we study graphically the behavior of the infectious population by reproducing the simulation showed in Fig. 1 of the author's article. We use NDSolve to generate a simulation of the behavior of both the infectious and susceptible population, and confirm the predictions for the new epidemic equilibrium point by looking at the final result of the numerical integration.

The Pulse Vaccination Strategy

Under the pulse vaccination scheme, a fraction of the population (hopefully much smaller than that of the constant vaccination strategy) is immunized, in single pulses (instead of in a constant and continuous way), between time intervals of $T$ years. This is reflected as an update to the rate of growth of the susceptible population, that will now decrease by a fraction of $p$ at each time of application of the pulsed vaccination, where $p$ is the proportion of vaccinated population. Mathematically, it is expressed as

(7)	$S(t_{n})=(1-p)S(t_{n-1})$,
	$t_{n+1} = t_{n} + T$.

Equation (7) should be added to the system of state equations describing the population in order to account for the periodic fluctuations in susceptible population.

The proportion of infectious population is assumed to be zero in order to explore an infection-free solution of the Pulse Vaccination SIR model ($I=0$, for all $t$).

Solving equation (7) for $S(t)$ and finding its stroboscopic map F such that $S_{n+1}=F(S_{n})$ (a map whose rule is to give the value of $S(t)$ after each pulse), it is possible to see that it has a single fixed point given by

(8) $S^{*}=F(S^{*})=\frac{(1-p)(e^{mT}-1)}{p-1+e^{mT}}$.

Noting that the orbit of map converges to this fixed point, the susceptible population $S(t)$ converges to a cycle of period T.

After solving equation (7) for $S(t)$ when $I(t) = 0, \forall{} t$, utilizing the fixed point value $S^{*}$, a complete infection-free solution is found:

(9)	$\tilde{S}(t) =$	$1 + \frac{pe^{mT}}{1-e^{mT}-p}e^{-m(t-t_{n})}, t_{n} \leq t \leq t_{n+1},$
		$S^{*}_{}, t = t_{n+1},$
	$\tilde{I}(t)=0$

The model may be linearized by setting the state variables as

$S(t)=\tilde{S}(t)+s(t),$ $I^{ }_{ }(t) = 0 + i(t)$,

where s and i are small perturbations of the population groups. Using this to expand Equation 1 in a Taylor series, the linearized equations result as (after neglecting higher-order terms)

(10) $\frac{ds}{dt} = -ms-\beta \tilde{S}(t)i,$ $\frac{di}{dt}=i(\beta \tilde{S}(t)-m-g),$

where $s(t_{n})$ is given by equation (7).

Equation (10) can be readily integrated over one pulse interval to yield

$i_{n+1}=i_{n}e^{\int^{t_{n+1}}_{t_{n}}(\beta \tilde{S}(t)-(m+g))dt}.$

Thus, $i_{n}$ will decrease exponentially fast if $\int^{t_{n+1}}_{t_{n}}(\beta \tilde{S}(t)-(m+g))dt < 0$.

This can be rearranged as

(11) $\frac{1}{T}\int^{T}_{0}\tilde{S}(t)dt<\frac{m+g}{\beta}=S_{c},$

where $S_c$ is a threshold for the mean value of $\tilde{S}(t)$ over a single pulse interval. If the value of $\tilde{S}(t)$ is below this threshold, the infection-free solution is stable. It can also be seen from equation 10 that as $i(t)$ approaches zero, $s(t)$ starts to converge to zero (damping the perturbations, making the infection-free equilibrium stable).

Reproducing Results

Using NDSolve, we numerically solve the system and then plot the results of simulating populational behavior as it converges to the expected equilibrium. We turn the rate equation for $S(t)$ into a piecewise function so as to account for the periodic perturbations caused by vaccination pulses. We divide the numerical integration procedure over several different time intervals spanning the period between each pulse, and join all pieces of the simulation together by use of Piecewise. The resulting integrated solution is then plotted over a time range spanning the entire interval of integration. This corresponds to the simulation shown in Figure 2 of the studied article.

Furthermore, we study the relation between the ratio of vaccinated population ($p$) and the maximum interpulse interval ($T_{\mathrm{max}}$). We use Solve to numerically solve the exact expression for the boundary for $T_{\mathrm{max}}$ under which the infection-free equilibrium is stable. We plot this numerical solution along with curves given by approximations through Taylor expansion and through an algebrical process described in . This corresponds to Figure 3 in the studied article.

Seasonally Enforced SIR model

Epidemiological data on measles suggest that the equilibria described so far are highly unlikely in practice. In fact, chaotic dynamics has been observed for measles epidemics in some cities. This demonstrates that the basic SIR model fails to capture fundamental aspects of measles dynamics that could account for more realism of the model's predictions.

A key feature of measles dynamics seems to be the seasonal variation in the infection rate. The infection rate might be affected by seasonal variation on the population's behavior, such as school period and immigration. The contact rate can then be made a function of time, varying along the year. This makes the system effectively three-dimensional, allowing for chaotic dynamics to arise.

The contact rate $\beta{}$ is given a base value $\beta{}_0$, and new parameter $\beta{}_1$ is created to account for the amplitude of seasonal variation. Under seasonal enforcing, the system of equations (1) should be added a new function to substitute the constant value of $\beta{}$, written below.

(11) $\beta{}(t)=\beta{}_0(1+\beta{}_1 cos 2\pi{}t), \text{ } \beta{}(t+1)=\beta{}(t)^{}_{}$

In fact, the system does display chaotic behavior if left to develop by itself. However, when a pulse vaccination is applied with useful parameters, these dynamics are dramatically altered and the system still converges to a periodic stable solution that leads to eradication of disease.

Updating the state equations we used to simulate the system under the pulse vaccination strategy, we then simulate the system accounting for variation in $\beta{}_0$, with different input values of $\beta{}_1$, the amplitude of variation. This is to demonstrate the capacity of seasonal variation to induce chaos and the pulse vaccination to suppress such chaotic behavior.

Mixed Vaccination Strategy

The mixed vaccination strategy consists of a rather straightforward implementation of both constant vaccination and pulse vaccination concepts in the SIR dynamical system. While it doesn't add too significantly to the richness of detail built so far, it nevertheless allows for a remarkable feature of the coupled constant and pulsed vaccinations. It significantly expands the maximum interpulse interval ($T_{\mathrm{max}}$), facilitating the efforts of holding the infectious population under a safe size limit that prevents epidemic outbreaks.

In order to demonstrate this feature, we again plot the ratio between the fraction of vaccinated population ($p$) and the maximum interpulse interval ($T_{\mathrm{max}}$), and show that the mixed vaccination in fact allows for much broader intervals between vaccination pulses. We use the same techniques we used to plot this relation for the pulse vaccination strategy alone.

Results

The Pure SIR Model

Stability Analysis

We first demonstrate the important relation of the constant $R_{0}$ with the stability of the system. The following mathematical development was performed in order to demonstrate it.

The jacobian matrix for equation (1) in the model description is

(12) $\left( \begin{array}{cc} -m-i \beta & -S \beta \\ i \beta & -g-m+S \beta \\ \end{array} \right)$

Taking one of the eigenvalues,

(13) $\lambda \equiv \frac{1}{2} \left(-g-2 m-i \beta +S \beta - \sqrt{\left(g^2-2 g i \beta -2 g S \beta +i^2 \beta ^2-2 i S \beta ^2+S^2 \beta ^2\right)}\right)$

substituting $I$ for $I^{*}$ and $S$ for $S^{*}$ and setting it to zero, we get

(14) $\lambda ^*\equiv \frac{1}{2} \left(-\frac{m \beta }{g+m}-\sqrt{m \left(4 (g+m)-4 \beta +\frac{m \beta ^2}{(g+m)^2}\right)}\right)=0.$.

If we suppose the square-root member is imaginary, by utilizing Reduce we find out that the condition for stability of the system is that one of the parameters be smaller than zero, what makes no biological sense. We proceed then on the assumption that the square-root member is real.

After algebraically treating (14), we get

$-\frac{m \beta }{g+m}=-\sqrt{m \left(4 (g+m)-4 \beta +\frac{m \beta ^2}{(g+m)^2}\right)}$

$\left(-\frac{m \beta }{g+m}\right)^2=m \left(4 (g+m)-4 \beta +\frac{m \beta ^2}{(g+m)^2}\right)$

$\frac{m \beta ^2}{(g+m)^2}=4 (g+m)-4 \beta +\frac{m \beta ^2}{(g+m)^2}$

$m \beta ^2=4 (g+m)^3-4 (g+m)^2 \beta +m \beta ^2_{ }$

$\beta =\frac{4 (g+m)^2}{m}+\frac{4 (g+m)^3}{m \beta }+\beta$

$4 (g+m)^2=\frac{4 (g+m)^3}{ \beta }$

$\beta =\frac{4 (g+m)^3}{4 (g+m)^2}$,

we arrive at the following condition,

$\beta =g+m^{ }_{ }.$

So that

(15) $R_0\left\{ \begin{array}{cc} >1 & \beta >g+m \\ <1 & \beta <g+m \\ \end{array} \begin{array}{c} \lambda ^* \text{ } \mathrm{is} \text{ } \mathrm{positive} \\ \lambda ^* \text{ } \mathrm{is} \text{ } \mathrm{negative} \\ \end{array} \begin{array}{c} \left(S^*,I^*\right) \mathrm{is } \text{ } \mathrm{unstable} \\ \left(S^*,I^*\right) \mathrm{is } \text{ } \mathrm{stable} \\ \end{array} \right.$

By observing the value of $R_0$ for the standard measles parameters (Table 3), we expect that, for this set of parameters, the trivial equilibrium point be unstable and the epidemic equilibrium point be stable.

Phase Plots

Making use of the techniques described in the past section, we generated plots to confirm our predictions on the latter paragraph. Figure 1 shows phase plots encompassing the entire meaningful range of the state variables. In Figure 2 we see that the nullclines seem to intercept each other at specifically $(S^*,I^*)$. This equilibrium point is highlighted on the left plot of Figure 2, while we demonstrate its stability by plotting a trajectory starting from a nearby point. The trajectory also helps getting a better feel for the flow of the system around the stable epidemic equilibrium point.

Figure 1. Global and local phase plot. Left: Global overview of the flow of the system. Both equilibrium points are pictured along what seem to be points where the nullclines intersect. Right: Local overview of the flow. A spiral flow towards the epidemic equilibrium point is hinted. [Used standard measles parameters]

Figure 2. Local view of the flow. Left: The intersection of the nullclines into the stable equilibrium is clearly shown. Right: A trajectory is plotted near the equilibrium point to clearly demonstrate its stability. [Used standard measles parameters]

The Constant Vaccination Strategy

Stability Analysis

We updated the rate of change of $S(t)$, set the rate equations to zero and using Reduce were able to find a new trivial equilibrium point $(1-p,0)$ as well as a relation for the new epidemic equilibrium point. Reduce returned that for the epidemic equilibrium point, $S(t) = S^*$ still. It gave us the following relation upon which we worked to find expressions for the new non-trivial equilibrium points.

	$I_c(t)=-\frac{m \left(\frac{g+m}{\beta }+p-1\right)}{g+m}$
	$I_c(t)=-\frac{m p}{g+m}-\frac{R_0 \left(m \left(\frac{g+m}{\beta }-1\right)\right)}{R_0 (g+m)}$
	$I_c(t)=-\frac{m p}{g+m}-\frac{m R_0 \left(\frac{1}{R_0}-1\right)}{\beta },\text{ } \text{as}\text{ } R_0=\frac{\beta }{g+m},\text{ } \text{and}\text{ } \text{finally}\text{ } \text{we}\text{ } \text{have}\text{ } \text{that}$
	$I_c(t)=\frac{m \left(1-R_0\right)}{\beta }-\frac{m p}{g+m}=I^*-\frac{m p}{g+m},\text{ } \text{given}\text{ } \text{that} \text{ }I^*=\frac{m \left(\frac{\beta }{g+m}-1\right)}{\beta }.$
	Thus, the equilibrium points under the constant vaccination strategy are
(16)	$(S^{*},I^{*} - \frac{mp}{g+m}), (1-p,0)$

The jacobian matrix for the SIR model under constant vaccination is still the same as (12), so, to assess stability, we plug in the new equilibrium points into the eigenvalues and compute the conditions for stability. Once more, if the square root member is imaginary, by using Reduce we find out that some parameters would have to take on negative values in order for the equilibrium point to be stable. We then solve the system on the assumption that the square-root member is real. After substituting the infectious equilibrium point values into the eigenvalue we get

$\left(\lambda _c\right){}^*\equiv \frac{1}{2} \left(\frac{\beta m (p-1)}{g+m}-\sqrt{m \left(\frac{\beta ^2 m (p-1)^2}{(g+m)^2}+4 (g+m)+4 \beta (p-1)\right)}\right)<0$

Since the first term must be negative (all parameters are positive and $p<1$) we need to invert the inequality after squaring. Then,

$\left(\frac{\beta m (p-1)}{g+m}\right)^2>\left(\sqrt{m \left(\frac{\beta ^2 m (p-1)^2}{(g+m)^2}+4 (g+m)+4 \beta (p-1)\right)}\right)^2$

$\frac{\beta ^2 m (p-1)^2}{(g+m)^2}>\frac{\beta ^2 m (p-1)^2}{(g+m)^2}+4 (g+m)+4 \beta (p-1)$

$0>-\beta +g+m+\beta p$

$p<\frac{\beta -g-m}{\beta }$

$p<1-\frac{1}{R_0}$

By use of Reduce we encounter the opposite relation for $p$ if we substitute the values of the trivial equilibrium. This confirms (4), (5) and (6) in the previous section. We also find that, for the standard measles parameters, $p_c$ is effectively 94.44%, what is really high.

Time Plot

Using the methods described in the previous section, we were able to reproduce Figure 1 of the original article (which corresponds to Figure 6 left), as well as plot its $S(t)$ counterpart. We confirm via numerical integration that the value of $S(t)$ at equilibrium remains unchanged, while that of $I(t)$ is reduced by the $- \frac{mp}{g+m}$ factor. The value of $I^*_c$ is $4.88769\times{}10^{-5}$, while that of the last time of integration is $4.86962.\times{}10^{-5}$.

Figure 3. Time series of both state variables. Both plots were generated during the same simulation. Until time $t=80\text{ }\mathrm{ years}$, the system is allowed to evolve unperturbed. However, at time $t=80\text{ }\mathrm{ years}$, a constant vaccination scheme is applied (equation (3)) with parameter $p=70%$ (fraction of vaccinated population). The system then converges to a new equilibrium. Left: Time series of $I(t)$. As predicted in (6), the normal equilibrium value of $I$ is reduced by a negative term. Right: Time series of $S(t)$. As predicted in (6), the normal value at equilibrium of $S$ doesn't change, and after a initial perturbation, $S$ settles again into the same equilibrium value $S^{*}$. [Used standard measles parameters]

The Pulse Vaccination Strategy

As described in the previous section, we divide the task of numerically integrating the rate equations into separate intervals spanning the value of the parameter $T$ (in years), use NDSolve to integrate the system of equations and join all the intervals in the plot with Piecewise. Using this technique, we produced Figure 4, which is a simulation of the system when it's close to equilibrium.

Instead of settling into a constant size, the susceptible population tends to grow with a determined amplitude until the next round of pulse vaccination takes away a fraction $p$ of susceptibles and turn them into recovered. The infectious population will also cycle, due to the intermittent pertubation of vaccination pulses. However, for certain values of $p$, the proportion of vaccinated population, that are very significantly smaller than that for constant vaccination, in combination with reasonably short $T$ periods between vaccination pulses, the size of the infectious population will drammatically decrease. This fact confirms the hypothesis that through pulsed vaccination policy, a smaller proportion of the population needs to be immunized.

Figure 4 shows two lines. The line red line corresponds to $S_c$, the threshold under which the average size of the susceptible population must remain during one interval in order for the system to remain on this stable state. The purple line shows $S^*$, the value $S(t)$ assumes immediately after each pulse of vaccination, and is a coordinate of the equilibrium.

Figure 4. Time series of $S$ and $I$ under the pulse vaccination strategy. We let the system stay at equilibrium until time $t=5$ and then start a pulse vaccination scheme with an interpulse interval $T=2\text{ }\mathrm{years}$. We used $p=51.6%$ (proportion of vaccinated population) rather than the 50% proposed in the original article only because it yielded results that were closer to the original figure. Left: Time series of $S(t)$ under the pulse vaccination strategy. $S_c$ represents the threshold for the average value of $S(t)$ in the cycle, above which the function $S(t)$ stops stably cycling. $S^*$ is the value of the function immediately after each pulse of vaccination over the stable range of the cycle. Right: Time series of $I(t)$ under the pulse vaccination strategy. It can be seen that the proportion of infectious population falls exponentially fast under the pulse vaccination strategy. In practical terms, we can say that it effectively reduces to zero, and so the epidemic disappears. [Used standard measles parameters]

Utilizing the criteria for stability of the infection-free equilibrium defined by (11), it is possible to derive an expression for the condition (an upper bound) for the value of $T$. This bound is called $T_{\mathrm{max}}$, and it represents the maximum length of the interval between pulses for which the infection-free equilibrium point is stable. This procedure yields the following expression,

(17) $\frac{m p T+\left(1-e^{m T}\right) (p-m T)}{m T \left(e^{m T}+p-1\right)}\leq \frac{g+m}{\beta }.$

It can only be solved numerically, for $T$ is expressed in terms of transcendental functions (functions for which it's not possible to solve for the independent variable, e.g. exponential functions).

Because of this, it's can be convenient to approximate this result. One way to do it is to make use of a Taylor expansion (this process consists of approximating the value of a function near a given point by taking sucessive derivatives of the function and using them to build a polynomial function that will approximate the original. Generally, terms containing higher derivatives are considered negligible, but the more terms a Taylor polynomial has, the better it will be as an approximation of the original function). After neglecting higher order terms, the process yields equation (18). The authors also derive another approximation in , which is represented by equation (19).

(18)	$T_{\mathrm{max}}\approx \frac{g p}{\beta m}\frac{1}{\left(-\frac{g}{\beta }-\frac{p}{2}+1\right)}$
(19)	$T_{\mathrm{max}}\approx \frac{1}{m}\log \left(\frac{p (g+m)}{\beta -g-m}+1\right)$

Figure 5 shows a plot of all three equations compared. These results corroborate the other part of the hypothesis, that if the time between vaccination pulses is smaller than some threshold, the epidemic will reach an acceptable stable state or disappear. Clearly, if a reasonable combination of the parameters $p$ and $T$ are used, the infection-free equilibrium point will be stable and the size of the infectious population will decrease exponentially fast as shown in Figure 4.

Figure 5. Blue represents equation (17), the exact result for the maximum value of $T$ for which the infection-free equilibrium is stable; Red represents equation (18), an approximation through Taylor expansions; Green represents equation (19), another approximation proposed in ; [Used standard measles parameters]

Seasonally Enforced SIR Model

When $\beta{}$ is allowed to vary, the systems ceases being two-dimensional and becomes three-dimensional, making possible for chaotic behavior to surface. In our simulations, we allowed the amplitude of seasonal variation, $\beta{}_1$, be adjusted freely. We first run the system while setting $\beta{}_1$ to zero (Figure 6).

Confirming that our simulation works for this established solution, we proceed to make $\beta{}_1=0.5$ (Figure 7). We observe chaotic patterns throughout our simulation of the state variables. Although the figure is static, very slight changes in the initial conditions can lead to very different outcomes at the final time and to different patterns of the curve.

As the infectious population grows unpredictably, risks that disease outbreaks occur are heightened. Notwithstanding, an infection-free equilibrium point still exists. It just takes to notice that the system (1) still has a solution where $I(t)=0$, and by looking at (1.1) we see that this solution makes $S(t)$ independent of seasonality (it makes the term containing $\beta{}(t)$ go to zero). Based on this we run our simulation one more time now applying pulses of vaccination for every two years (Figure 8). We see that this equilibrium will be attained, at that pulse vaccination changes the dynamics of seasonally enforced SIR drastically.

Figure 6. Simulation of the system with amplitude of seasonal variation parameter $\beta{}_1$ set to zero. The system effectively behaves as an unperturbed/unaffected pure SIR system. Please note that these plots are analogous to those of Figure 3. [Used standard measles parameters]

Figure 7. Simulation of the system with amplitude of seasonal variation $\beta{}_1$ set to 0.5. Middle: $\beta{}(t)$ now varies covering a certain amplitude each year. Left and Right: Clearly making $\beta{}$ a function of time has induced the system to behave chaotically. [Used standard measles parameters]

Figure 8. Simulation of the system with amplitude of seasonal variation $\beta{}_1$ set to 0.5 and under pulse vaccination with interpulse interval $T=2\text{ }\mathrm{years}$ and vaccinated population proportion $p=51.6%$. Chaotic behavior is suppressed and the system resembles much that of Figure 4. [Used standard measles parameters]

Mixed Vaccination Strategy

In the mixed vaccination strategy, the same equilibrium of the pulse vaccination scheme coexist, except that, as it was for the constant vaccination strategy alone, the value in the infection-free equilibrium of $S(t)$ will be reduced by a factor of $(1-p_c)$, where $p_c$ is the proportion of population that will be constantly vaccinated. This means that $S^*$ for the pulse vaccination strategy will be multiplied by $(1-p_c)$. Since $S^*$ was used as a condition for stability (recall its formula from Table 3 and observe equation (17)), the new exact solution for $T_{\mathrm{max}}$ is just (17) multiplied by $(1-p_c)$ (equation (20)). In Figure 9 we use the same techniques we used to produce Figure 5 to plot this relation between $p$ for the pulse vaccination and $T_{\mathrm{max}}$, keeping $p_c$ constant at 85%. We also plot an approximation proposed in (equation (21)) as well as the exact solution for the pulse vaccination strategy alone.

(20) $\frac{\left(m p T+\left(1-e^{m T}\right) (p-m T)\right)\left(1-p_c\right)}{m T \left(e^{m T}+p-1\right)}<\frac{g+m}{\beta }$

(21) $T_{\max }= \frac{p g}{\beta m \left(\left(1-p_c\right)(1-p/2)-g/\beta \right)}$

Clearly, if feasible, the mixed vaccination scheme turns out to be the most flexible and safest of the strategies studied so far.

Figure 9. Red represents equation (20), the exact result for the maximum value of $T$ for which the infection-free equilibrium is stable under the mixed vaccination policy with $p_c = 85%$; Green represents equation (21), an approximation proposed in ; Blue represents equation (17), the exact result for the maximum value of $T$ for which the infection-free equilibrium is stable under only the pulse vaccination scheme. [Used standard measles parameters]

Discussion

In our reevaluation of the model, we were able to reproduce some algebraic results and the qualitative results rather faithfully. Given the theoretical character of the study, the mechanistic approach established by the SIR model, the reasonable parameters the authors utilize and the assumptions and conditions under which they develop their analysis, our results demonstrate that the initial expectations proposed in the introduction should be satisfactorily fulfilled.

The hypotheses centered about the efficacy of pulse vaccination as a strategy to combat epidemics, particularly as compared to constant vaccination. They state that a smaller proportion of the population needs to be immunized under the pulse vaccination strategy and that as long as the interval between pulses doesn't exceed a limit, the epidemic will disappear. The first statement was confirmed by us. While in the pure/unaffected SIR model the system tends to a stable epidemic equilibrium, under pulse vaccination the fraction of susceptible population floats within a stabe cycle if its size averages under a certain limit. This is achieved by significantly smaller proportions of vaccinated population if compared to the proportions needed for constant vaccination alone.

The other part of the hypothesis was also confirmed through plotting and unveils one weakness of this study. In figure 5, we plot the dependence of the maximum interpulse interval on the proportion of vaccinated population. In simple terms, if we use any value for the interpulse interval below those curves, the epidemic will be extinguished. This was portrayed in figure 4, where the proportion of infectious population decays exponentially fast. The apparent flaw in the model is that it doesn't account for the discreteness of the system, it rather treats, in all of its simplicity, its components as elements of continuous nature. For instance, if we take a population of million, the SIR model will still compute values for the infectious population smaller (much smaller) than $10^{-6}$. It fails to consider that instead of decreasing to still smaller values, the infectious population should effectively be reduced to zero, and the epidemic should never resurface unless some immigration of infected individuals occur. This is due, we're inclined to think, to a mere intention of making an exploratory analysis of the system that might give clues on how to extend, add complexity and realistic attributes to mechanistic models of disease dynamics in order to improve their descriptive and predictive capabilities. The model also fails to depict other features of disease dynamics, such as the inherent heterogeneity of the interactions between susceptible and infected populations (the SIR model generally assumes mass-action dynamics to account for rather homogeneous interactions among populations). Nevertheless, even though being one of the simplest, the SIR model is among the most accurate of all biological models.

This is a developing field, and the authors themselves explore such additions to the model in order to explore its realistic capabilities. For instance, in an earlier paper , they do the same description of the model but add more elements to the analysis, such as seasonal forcing and mixed strategies (concomitant pulsed and constant vaccination procedures). These additions were the extensions we made on this paper. We successfully demonstrated that the SIR model can account for seasonal variation and chaotic dynamics. We also used this to strengthen our confidence in the efficiency of vaccination strategies.

The SIR model is rather a staple on the study of disease dynamics, and much has been added and built over its essential propositions. Lu et al. extend over the model destroying one of its assumptions that is that there is no vertical transmission (parents don't contaminate offspring through birth or gestation). In that work, Lu et al. also go on to evaluate the theoretical effectiveness of the different strategies and add to the discussion that one must take into the account the different nature of the vaccines, implying that there may be circunstances where a constant vaccination strategy would be more reasonably appliable. Clearly, there's much to develop over the rather humble foundations of the basic SIR model.

All in all, we found no discrepancies between our replication and the original article itself. In spite of its numerous but obvious limitations, such as being unconsiderate to discreteness of data, unrespectful to alternative processes of disease transmission and even ignorant to seasonal variation of transmission rates, it still forms a strong base over which to analyze and plan social, political, economical and scientific efforts to combat epidemics.

Appendix

My Mathematica notebook

Code for Seasonal Enforcing Simulations