The Kermack-McKendrick Model of Epidemics and Its Extensions

Jan Hurt and Daniel Hlubinka

Charles University of Prague, Czech Republic

1 Introduction

The classical deterministic model of epidemic was proposed by Kermack and McKendrick in 1927, see [4]. The model is considered to be appropriate for highly infectious diseases with a fast recovery and/or fatality. Influenza, measles, chickenpox, mumps, and scarlet fever are examples of such diseases. It assumes a constant rate of infection and a constant recovery rate. There are various extensions of the model in that the rate of infection is not constant, the size of population varies in time etc. Here we mention two practically important cases; intervention either in the form of government measures or vaccination and a random size of population.

2 Deterministic models

2.1 Classical model

It is assumed that the population of size N, constant in time t, consists of three parts: (Subscript[X, t]), the number of susceptibles (exposed to the infection), (Subscript[Y, t]), the number of infectives, and (Subscript[Z, t]), the number of removals (recovered, dead, or isolated), (Subscript[X, t]) + (Subscript[Y, t]) + (Subscript[Z, t]) = N. The model is governed by the system of ODE (the derivatives are taken with respect to t)

\(X\& . \_t = \(-βX\_t\) Y\_t, \ 
    X\_0 = x\_0, 
Y\& . \_t = \(βX\_t\) 
          Y\_t - γY\_t, 
    Y\_0 = y\_0, 
Z\& . \_t = γY\_t, \ 
    Z\_0 = \(\(z\_0\)\(.\)\(\ \ \ \ \)\)\)

This model is also known as an SIR model (Susceptibles, Infected, Removed). Positive constants β and γ represent the infection and recovery rates, respectively. Thus the transition X⟶Y is realized with rate β, Y⟶Z with rate γ, consequently X⟶Z with rate ρ := γ/β, the so called relative removal rate. The dimension of β is (((time × size of population))^(-1)), that of γ (((time))^(-1)), and that of ρ is size of population. In what follows we will measure time in days. Without any loss of generality, the size of population will be set to one: N:=1.

For β > 0, γ > 0, (Subscript[x, 0]) + (Subscript[y, 0]) = N, (Subscript[x, 0])>0 , and (Subscript[y, 0])>0 there exists a solution to the above system. The solution ((Subscript[X, t]), (Subscript[Y, t]), (Subscript[Z, t])) is nonnegative, the number of susceptibles (Subscript[X, t]) is nonincreasing, the number of removals (Subscript[Z, t]) is nondecreasing. As t⟶∞ there exist positive limits (Subscript[X, ∞]) and (Subscript[Z, ∞]) of (Subscript[X, t]) and (Subscript[Z, t]), respectively. If ρ<(Subscript[x, 0]) then there exists the global maximum of infectives, (SubPlus[Y])=max (Subscript[Y, t]), say.

2.1.1 Example

To get an idea about the magnitudes of the parameters involved, we give the expert's guesses for a (possibly pandemic) influenza. The choice β=0.6/N means that the probability of the transfer of the disease from the infected individual to a susceptible is approximately 0.6. The choice γ = 0.2 means that approximately 1/γ = 5 days are needed for the transition from infected to removals.

2.1.1 The code

Since the given system of ODE may only be solved numerically, we make use of Mathematica NDSolve procedure. The time interval of interest is [t1, t2], n0 is the initial size of population, the other symbols have been explained above. For the future use we also introduce variable z0, the initial number of removed. The output is the list, the first part of which is the plot of the solution and the second part is the list of interpolating functions x, y, z.

Clear[SIR];
Needs["Graphics`"];
SIR[β_, γ_, x0_, z0_, n0_, {t1_, t2_}, 
        plotlabel_: "Evolution"]:=
Module[{solsys, x, y, z},
solsys = 
              NDSolve[{x'[t] == -β x[t] y[t], y'[
                      t] == β x[t] y[t] - γ y[
                          t], z'[t] == γ y[t], 
                  x[t1] == x0, y[t1] == n0 - (x0 + z0), 
                  z[t1] == z0}, {x, y, z}, {t, t1, t2}];
{x, y, z} = {x, y, z} /. 
                Flatten[solsys];
{Plot[{x[t], y[t], z[t]} // Evaluate,{t, t1, t2},PlotRange -> All,
PlotStyle -> {Dashing[{0.005, 0.01}], GrayLevel[0], Dashing[{0.03, 0.03}]},
\(TextStyle -> {FontFamily -> "Tahoma", FontSize -> 12}\),
PlotLabel -> plotlabel,
AxesLabel -> {" t (in \
days)", " "},
\(PlotLegend -> {"susceptibles", "infected", "removed"}\),
LegendPosition -> {0.8, 
                        0.1}],{x, y, z}}]

This illustrates the solution with the data from Example 2.1.1:

example211 = 
        SIR[0.6, 0.2, 0.999, 0, 1, {0, t∞ = 75}];

The following plot shows the infectives only. We see that some 30 per cent of infectives in the worst case means the pandemic.

y211 = example211[[2]][[2]]  (* = y[
              t]  ...  number of infectives *) 
Plot[y211[t], {t, 0, t∞}];

2.2 Models with changing infection rate

So far we assumed that the infection rate β is a constant. Nevertheless, there are many situations appearing in practise where such an example is unrealistic. The infection rate may be considered either as a function of time or of number of removals. We mention two particular cases.

2.2.1 Kendall's model

In this model it is supposed that the epidemic is governed by the classical model but the infection rate changes with the number of the removals. In other words, β is a function of z. The so called  Kendall's parameter β[z] which is given by
\(β[
          z] = \(2  β\_0\)\/\(\((1 - \(z\ β\_0\)\/γ)\) + \((1 - \(z\ β\_0\)\/γ)\)\^\(-1\)\), \ 
      0 ≤ z < \(\(γ\/β\_0\)\(.\)\)\)
Note that β[z] is decreasing as the number of removals is increasing and hence it may be used to model epidemics which are loosing the strength with decreasing number of susceptibles. Another important fact is that there exists an explicit solution to the classical deterministic SIR model with Kendall's β[z]. In our case, however, we use this function in Section 3.

2.2.2 Intervention models

The intervention usually takes place if the infection spreads to an unacceptable level.The intervention may be incorporated into the model by changing the parameter β. Since there is almost always to solve the equations numerically, we can model the intervention by any reasonable function. Moreover, by Mathematica we are also able to model some strategic decisions concerning spreading the disease. 

In principle, the intervention may take two forms or both of them. The first consists in restrictions in contacts (closing schools, cinemas, etc.), the second one in vaccination of a part of a population (usually of people with intensive contacts with other people like drivers of the public transport, doctors, shop assistants) or just the whole population. Vaccination of a part of population would eventually lead to a reduction of susceptibles. All these measures result in changing β. 

Here we will develop a model in which the intervention will take place if the number of infectives reaches a prescribed percentage of population p, say, at time (Subscript[t, p]). From the technical point of view, at time (Subscript[t, p]) we decrease β. Since under our assumptions y[t] possesses just one global maximum on [0, ∞], we will seek the first passage time of y[t] through Np. This and other actions described later can be surely done by very powerful method EventLocator for NDSolve but the authors of this contribution were not able to realize it in a reasonable time. We attack the problem in the following way. First we solve the problem on a sufficiently long time interval. Second we locate the first time where the maximum is achieved and then we are looking for (Subscript[t, p]). Such an approach almost automate the procedure.

We will illustrate the procedure on the Example 2.1.1 step by step until we present the general procedure.

2.2.2.1 Step by step

Remember that y211 is the interpolating function of infectives. First we find the approximate maximum and later use it as the initial solution to the differential equation y'[t]=0 which gives the unique maximum.

yaux211 = 
      Max[Table[
          y211[t], {t, Floor[t∞]}]]
taux211 = Position[Table[y211[t], {t, Floor[t∞]}], 
          yaux211][[1, 
        1]]
ymax211 = 
      FindRoot[y211'[t] == 0, {t, 
          taux211}]
y211[t] /. ymax211

Next we prescribe the upper percentage level of the infectives p and we get the desired starting time of the intervention. Due to monotonicity we choose the initial solution 90 per cent of the maximum.

p = 0.1;
tp = 
      t /. FindRoot[y211[t] == p, {t, 0.9  t /. ymax211}]

We  can look what will happen if we change β → 0.4 at this time. First we must extract the final values of x, y, z at (Subscript[t, p]).

{xp, yp, 
        zp} = {example211[[2]][[1]][
          tp], example211[[2]][[2]][
          tp], example211[[2]][[3]][tp]}

example211before = 
        SIR[0.6, 0.2, 0.999, 0, 
          1, {0, tp}, "Evolution before intervention"];
example211after = 
        SIR[0.4, 0.2, xp, zp, 
          1, {tp, t∞}, "Evolution after intervention"];

The following plot shows the infectives only. We see that now the maximum of the infectives was reduced to some 17 per cent.

Show[
        Plot[example211before[[2]][[2]][t], {t, 10, tp}, 
          DisplayFunction -> Identity], 
        Plot[example211after[[2]][[2]][t], {t, tp, 
            25 (*t∞*) }, DisplayFunction -> Identity], 
        DisplayFunction -> $DisplayFunction];

2.2.2.2 General procedure

Clear[SIRintervention];
Needs["Graphics`"];
\(SIRintervention[p_, β0_, β1_, γ_, 
        x0_, z0_, 
        n0_, {t1_, t2_}] := 
Module[{outSIR, outSIR1, 
          outSIR2, yaux211, taux211, ymax211, y211, solsys, y211before, 
          y211after, t∞}, 
t∞ = 
          t2; 
y211 = \(\((solsys = 
                  SIR[β0, γ, x0, z0, 
                    n0, {t1, 
                      t∞}])\)[[2]]\)[[2]]; 
yaux211 = Max[Table[y211[t], {t, Floor[t∞]}]]; 

        taux211 = \(Position[Table[y211[t], {t, Floor[t∞]}], 
              yaux211]\)[[1, 1]]; 

        ymax211 = 
          FindRoot[\(y211'\)[t] == 0, {t, 
              taux211}]; 
tp = 
          t /. FindRoot[
              y211[t] == p, {t, 
                0.9  t /. ymax211}]; 
{xp, yp, 
            zp} = {\(\(solsys[[2]]\)[[1]]\)[t] /. 
              t -> 
                tp, \(\(solsys[[2]]\)[[2]]\)[t] /. 
              t -> 
                tp, \(\(solsys[[2]]\)[[3]]\)[t] /. 
              t -> tp}; 
y211before = 
          SIR[β0, γ, x0, z0, 
            1, {0, tp}, "Evolution before intervention"]; 
y211after = 
          SIR[β1, γ, xp, zp, 
            1, {tp, t∞}, "Evolution after intervention"]; 
Show[
          Plot[\(\(y211before[[2]]\)[[2]]\)[t], {t, 0.5\ tp, tp}, 
            DisplayFunction -> Identity], 
          Plot[\(\(y211after[[2]]\)[[2]]\)[t], {t, tp, 
              2.5\ tp (*t∞*) }, DisplayFunction -> Identity], 
          DisplayFunction -> $DisplayFunction];
		]\)

SIRintervention[0.1, 0.6, 0.4, 0.2, 0.999, 0, 
      1, {0, t∞ = 100}]

2.3 Model with vaccinated subpopulation

The vaccination is a powerful tool in fighting infective disease. There are two possibilities, which can also be combined.

There is some percentage (Subscript[p, v]) of people vaccinated before the epidemics starts and are therefore immune. It is the same model as before with the exception that now (Subscript[z, 0  v]) = (Subscript[p, v]) (Subscript[n, 0]) and (Subscript[x, 0  v]) = (Subscript[x, 0])(1 - (Subscript[p, v])).

The vaccination is applied after the epidemics starts and there is some vaccination rate v which describes the speed of transition between the susceptibles and removals without passing the infected stage. The vaccination rate can be changing depending on the number of infectives or removals. Hence it is quite different model.

2.3.1 Vaccination with rate v

A new subpopulation is considered in this model. Beside the susceptibles, infectives and removals the subpopulation of  vaccinated, denoted  W, appears. The system of the differential equations is changed to

\(X\& . \_t\)=-β(Subscript[X, t]) (Subscript[Y, t])-(Subscript[vX, t]),          (Subscript[X, 0]) = (Subscript[x, 0]),
\(Y\& . \_t\)=β(Subscript[X, t]) (Subscript[Y, t])-(γ + v)(Subscript[Y, t]),(Subscript[Y, 0]) = (Subscript[y, 0]),
\(Z\& . \_t\)=(γ + v)(Subscript[Y, t]),                 (Subscript[Z, 0]) = (Subscript[z, 0]),
(Subscript[\(W\& . \), t])= (Subscript[vX, t]),                             (Subscript[W, 0]) = (Subscript[w, 0]).

Description of the model:
the rate v is the proportion of the population which can be vaccinated every day. Both infectives and susceptibles are subject to the vaccination but while susceptibles become vaccinated (denoted W  here), the infectives become removed — they actually have passed the disease.

2.3.2 Example with no vaccinated individuals at t = 0

Clear[x, y, z, w];
n0 = 1; 
    x0init = 0.999; β0 = 0.6; γ0 = 0.2; t∞ = 240; 
    v0 = 0.005;
solsys = 
        NDSolve[{x'[t] == -β x[t] y[t] - 
                  v x[t], y'[
                  t] == β x[t] y[t] - γ y[t] - 
                  v y[t], z'[t] == γ y[t] + 
                  v y[t], w'[t] == v x[t], x[0] == x0, 
              y[0] == n0 - x[0], z[0] == 0,  
              w[0] == 0} //. {n -> n0, 
              x0 -> 
                x0init, β -> β0, γ -> γ0, 
              v -> v0}, {x, y, z, w}, {t, 0, 
            t∞}];
x = (x /. 
            solsys)[[1]]; 
    y = (y /. solsys)[[1]]; 
    z = (z /. solsys)[[1]]; 
    w = (w /. 
            solsys)[[1]];

Plot[{x[t], y[t], z[t], w[t]} // Evaluate,{t, 0, 90},PlotRange ->  All,PlotStyle -> {Dashing[{ .005, 0.01}], GrayLevel[0], 
              Dashing[{ .03, 0.03}], Dashing[{0.03, 0.01, 0.0025, 0.01}]},TextStyle -> {FontFamily -> "Tahoma", 
              FontSize -> 12},
PlotLabel -> "Vaccination I",
AxesLabel -> {" t (in days)", " \
"},
\(PlotLegend -> {"susceptibles", "infected", "removed", "vaccinated"}\),
LegendPosition -> {0.8, 0.1}];

A vaccination of a slow rate v applied after the beginning of the epidemic cannot help very much. Only very fast vaccination can have some important effect in such a case.

The combination of the two type of vaccination is, however, quite powerful tool.  The results of the preventive vaccination are fascinating; With only 15 per cent of the population vaccinated, the total number of infected people is decreased substantially. The direct vaccination has almost no effect unless the vaccination rate is high (and therefore a relatively big proportion of the population is vaccinated every day).

2.3.3 Example with initially vaccinated individuals at t = 0

Consider the preventive vaccination of 15 percent of the population, (Subscript[p, v]) = 0.15. The vaccination rate after the epidemic starts is relatively small, v = (3/1000); every day we vaccinate three out of one thousand of susceptibles and infectives.

Clear[x, y, z, w];
n0 = 1; pv = 0.15;  
    x0init = (1 - pv) 0.999; 
    w0init = pv n0;  β0 = 0.6; γ0 = 0.2; t∞ = 200; 
    v0 = 0.003;
solsys = 
        NDSolve[{x'[t] == -β x[t] y[t] - 
                  v x[t], y'[
                  t] == β x[t] y[t] - γ y[t] - 
                  v y[t], z'[t] == γ y[t] + 
                  v y[t], w'[t] == v x[t], x[0] == x0, 
              w[0] == w0, y[0] == n0 - x[0] - w[0],  
              z[0] == 0} //. {n -> n0, x0 -> x0init, 
              w0 -> 
                w0init, β -> β0, γ -> γ0, 
              v -> v0}, {x, y, z, w}, {t, 0, 
            t∞}];
x = (x /. 
            solsys)[[1]]; 
    y = (y /. solsys)[[1]]; 
    z = (z /. solsys)[[1]]; 
    w = (w /. 
            solsys)[[1]];

Plot[{x[t], y[t], z[t], w[t]} // Evaluate,{t, 0, 90},PlotRange ->  All,PlotStyle -> {Dashing[{ .005, 0.01}], GrayLevel[0], 
              Dashing[{ .03, 0.03}], Dashing[{0.03, 0.01, 0.0025, 0.01}]},
\(TextStyle -> {FontFamily -> "Tahoma", FontSize -> 12}\),
PlotLabel -> "Preventive \
Vaccination ",
\(AxesLabel -> {" t (in days)", " "}\),
\(PlotLegend -> {"susceptibles", "infected", "removed", "vaccinated"}\),
LegendPosition -> {0.8, 0.1}];

3 A stochastic model

This part will be devoted to a stochastic generalization in that the population size is a stochastic process taking into account immigration and/or emigration, e. g. The population size (Subscript[N, t]) = (Subscript[X, t]) + (Subscript[Y, t]) + (Subscript[Z, t]) + (Subscript[W, t]) is governed by the stochastic differential equation (also called  Engelbert-Schmidt equation)

d(Subscript[N, t]) = σ[(Subscript[N, t])] (Subscript[N, t]) d(Subscript[B, t]),  
    (Subscript[N, 0]) = (Subscript[n, 0]),

where σ[ . ] is a bounded function with a compact support and (Subscript[B, t]) is a standard Brownian motion - Wiener process. The processes (Subscript[X, t]), (Subscript[Y, t]), (Subscript[Z, t]), and (Subscript[W, t]) behave like in the deterministic case, in average. Various approaches will be analyzed and a simulation model will be developed.

The infection parameter β will now be considered as a function of the removals (Subscript[Z, t]) instead of being constant. The SIR model with vaccination and stochastic population size can be written as the system of  four stochastic differential equations

d(Subscript[X, t])=-β[(Subscript[Z, t])](Subscript[X, t])(Subscript[Y, t])dt-v (Subscript[X, t])dt+(Subscript[X, t])σ[(Subscript[N, t])]d(Subscript[B, t]),           (Subscript[X, 0]) = (Subscript[x, 0]),
d(Subscript[Y, t])=β[(Subscript[Z, t])](Subscript[X, t])(Subscript[Y, t])dt-(γ + v) (Subscript[Y, t]) dt+(Subscript[Y, t]) σ[(Subscript[N, t])] d(Subscript[B, t]),  (Subscript[Y, 0]) = 
        (Subscript[y, 0,])
\(dZ\_t = \(\((γ + v)\) \(Y\_t\) dt + \(Z\_t\) σ[
              N\_t] dB\_\(\(t\)\(,\)\(\ \)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Z\_0 = 0\), 
dW\_t = 
      v\ \(X\_t\) dt + \(W\_t\) σ[
            N\_t] dB\_t, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ W\_0 = 
      w\_\(\(0\)\(,\)\)\)

where (Subscript[N, t]) = (Subscript[X, t]) + (Subscript[Y, t]) + (Subscript[Z, t]) + (Subscript[W, t]), γ and v have the same meaning as above.

Apparently, the numerical solution to the system requires to simulate four diffusion processes (Subscript[X, t]),  (Subscript[Y, t]),  (Subscript[Z, t]), and (Subscript[W, t]), respectively. In fact it is sufficient to simulate just one diffusion process (Subscript[N, t]), since the stochastic differential equations above may be rewritten as a system of ODE with random coefficients

dx (t)=-α[z (t), (Subscript[N, t])] x (t) 
              y (t) dt-v (Subscript[X, t])dt,                   x (0) = ((Subscript[x, 0])/(Subscript[N, 0])),
dy (t) = α[z (t), (Subscript[N, t])] 
              x (t) 
              y (t) dt - (γ + 
                  v) y (t) dt,   y (0) = ((Subscript[y, 0])/(Subscript[N, 0])),
dz (t)  =  (γ + 
                v) y (t) dt,                                                               z (0) = 0,
dw (t)=v (Subscript[X, t])dt,                                                                                    w (0) = ((Subscript[w, 0])/(Subscript[N, 0])),

where α(z, n)=n β(n·z), and (Subscript[N, t]) is given by the Engelbert-Schmidt equation. In particular case, it is possible to use the deterministic solution of the ODE described in Section 2 and multiply the solution ((Subscript[X, t]),  (Subscript[Y, t]),  (Subscript[Z, t]),  (Subscript[W, t])) by the simulated population size (Subscript[N, t]). Such a special case simplifies calculations and analysis a bit. Usually we have to solve a system of differential equations with random coefficients (in particular the coefficient  α depends on the population size) numerically. In both cases, the first step is to simulate the population size.

3.1 The population size

3.1.1 Step I: Engelbert-Schmidt equation of the population size

The strong solution of the Engelbert-Schmidt equation

d(Subscript[N, t]) = σ[(Subscript[N, t])] (Subscript[N, t]) d(Subscript[W, t]),  
    (Subscript[N, 0]) = (Subscript[n, 0])

exists if σ is a Lipschitz function. We typically choose

σ[n]={({{ \(\(\(\((n - 
                                        a)\)\^ν\_1\) \((b - n)\)\^ν\_2\)\/\(\(n\^ν\_3\) \@\(a + b\)\)\)        if a ≤n≤b,}, {0   			  otherwise,}})

where a < b are the bounds of the population size, and  σ in this case need not to be locally Lipschitz. However, there still exists a weak solution to the Engelbert-Schmidt equation if

\(∫\_\(-ε\)\%ε\(\((σ[
                  x + n])\)\^\(-2\)\) dx = \(∞\ \ \ \(∀ ε > 
              0\ \ \ ⟺\ \ \ σ[n]\) = 
        0. \)\)

The weak solution can be constructed in a following way. Let  {(Subscript[B, t]),  t ≥ 0} be a standard Brownian motion defined on some complete probability space. Denote

(Subscript[A, t])=((∫_0)^t)((((Subscript[B, u])·σ[(Subscript[B, u])]))^2)du

a strictly nondecreasing process, and

(Subscript[τ, t]) = inf {s : (Subscript[A, s]) > t}

the time of the first entry of (Subscript[A, s])  to (t, ∞). According to Engelbert-Schmidt, the process

(Subscript[N, t]) = (Subscript[B, (Subscript[τ, t])])

is the unique weak solution to the SDE for the population size.

3.1.2 Step II: Simulation of the population size

We start with simulated standard Brownian motion (Subscript[B, s]) up to its first entry to a or b. Let us consider a sequence of independent N(0,1) distributed r.v. (Subscript[Z, i]) and define simulated Brownian motion with reference points

\(0 = \(t\_0 < 
          t\_1 < ⋯ < τ\_{a, b} = 
        inf {\(t\_j\) : \ B\_\(t\_j\) ≤ a, \ 
            or\ \ B\_\(t\_j\) ≥ 
              b}\), 
B\_\(t\_i\) = 
      B\_\(t\_\(i - 1\)\) + \(\@\(t\_i - t\_\(i - 1\)\)\) Z\_i, \ 
    B\_0 = n\_0, 
B\_u = B\_\(t\_\(i - 1\)\) + \(\(u - t\_\(i - 1\)\)\/\(t\_i - 
                t\_\(i - 1\)\)\) \((B\_\(t\_i\) - B\_\(t\_\(i - 1\)\))\), \ 
    t\_\(i - 1\) ≤ u ≤ t\_i, \)

where the time step for simulations of a standard Brownian motion is usually considered fixed, Δ := (Subscript[t, i]) - (Subscript[t, i - 1]). Using this simulated piecewise linear process we calculate the increasing process (Subscript[A, t]) as piecewise linear with

(Subscript[A, (Subscript[t, i])])=(Subscript[A, (Subscript[t, i - 1])])+(((Subscript[t, i]) - (Subscript[t, i - 1]))/(Subscript[Z, i]))(∫_((Subscript[B, (Subscript[t, i - 1])]))^((Subscript[B, (Subscript[t, i])])))(((x·σ[x]))^(-2))dx, (Subscript[A, 0]) = 0

or more simply by Riemann sums

\(A\_\(t\_i\) = 
      A\_\(t\_\(i - 1\)\) + \(t\_i - t\_\(i - 1\)\)\/\((B\_\(t\_\(i - 1\)\)·σ[B\_\(t\_\(i - 1\)\)])\)\^2, \ 
    A\_0 = 0. \)

Then we use (Subscript[N, (Subscript[A, (Subscript[t, i])])]) = (Subscript[B, (Subscript[t, i])]) and we employ some interpolation (may be piecewise linear) to calculate (Subscript[N, u]) for any u.

3.1.2.1 Remark

A simple simulation of the Brownian motion using  (Subscript[t, i]) - (Subscript[t, i - 1]) = Δ has a major disadvantage. The increments of  (Subscript[A, (Subscript[t, i])]) - (Subscript[A, (Subscript[t, i - 1])]) become extremely large as B approaches the limit a or b, and on the other hand they can be very small if B is far from a and b. Hence the simulated time differences for the population size process may differ substantially.
	Some results from the classical theory like Lévy modulus of continuity or Law of iterated logarithm can be used to construct self-adjusted time increments,  (Subscript[t, i]) - (Subscript[t, i - 1]) = ξ[(Subscript[B, (Subscript[t, i - 1])])], for the simulated Brownian motion.

For the sake of simplicity we prefer to adjust the time increment of the simulated Brownian motion with the aim to control the increments of (Subscript[A, t]).  Hence the simulated time step of (Subscript[N, t]) is kept in reasonable bounds. It  For the numerical solution of the epidemic model it is quite useful since we need to interpolate the simulated process.

3.1.3 Example

We consider the initial population size (Subscript[N, 0]) = 1025, the bounds a = 1000,  b = 1050, and the diffusion parameter σ in the form

σ[n]={({{3 \(\(\(\((n - a)\)\^\(1/
                                      2\)\) \((b - n)\)\^\(1/2\)\)\/\((a + b)\)\^\(1/2\)\)        if a ≤n≤b,}, {0   			  otherwise.}})

Initial values:

Needs["Statistics`"]
SeedRandom[13];
aA = 1000; bB = 1050; n0 = ((aA + bB)/2); 
    n = 1000; δ = (1/n);
distr = 
        NormalDistribution[0, (Sqrt[δ])];
sigma2star[
        x_] := ((( 3 ((Sqrt[(x - aA) (bB - x)])/(((aA + bB))^(1/2)))))^2)

3.1.3.1 Simulation of the Wiener proces

SeedRandom[13];
n = 10; δ = (1/n); 
    distr = NormalDistribution[0, (Sqrt[δ])];
walk11 = 
        NestWhileList[(# + Random[distr] &), n0, 
          aA + 5 < # < bB - 5 &];
times1 = 
        FoldList[Plus, 0, 
          Table[(1/n), {Length[walk11] - 1}]];
n = 100; δ = (1/n); 
    distr = NormalDistribution[0, (Sqrt[δ])];
walk12 = 
        NestWhileList[(# + Random[distr] &), Last[walk11], 
            aA + 2 < # < bB - 2 &] // Rest;
times2 = 
        FoldList[Plus, Last[times1], Table[(1/n), {Length[walk12]}]] // 
          Rest;
n = 1000; δ = (1/n); 
    distr = NormalDistribution[0, (Sqrt[δ])];
walk13 = 
        NestWhileList[(# + Random[distr] &), Last[walk12], 
            aA + 0.1 < # < bB - 0.1 &] // Rest;
times3 = 
        FoldList[Plus, Last[times2], Table[(1/n), {Length[walk13]}]] // 
          Rest;
n = 10000; δ = (1/n); 
    distr = NormalDistribution[0, (Sqrt[δ])];
walk14 = 
        NestWhileList[(# + Random[distr] &), Last[walk13], 
            aA < # < bB &] // Rest;
times4 = 
        FoldList[Plus, Last[times3], Table[(1/n), {Length[walk14]}]] // 
          Rest;

globalwalk = 
        Join[walk11, walk12, walk13, walk14];
globaltimes = 
        Join[times1, times2, times3, times4];
Map[Length, {globaltimes, globalwalk}];
wienerwalk = 
        Transpose[{globaltimes, globalwalk}];
\(\(ListPlot[wienerwalk, 
        PlotJoined -> 
          True, 
TextStyle -> {FontFamily -> "Tahoma", 
            FontSize -> 
              14}, 
PlotLabel -> "Wiener process"];\)\)

3.1.3.2 Calculation of the process (Subscript[A, t])

timesdifferences = 
        Rest[globaltimes] - 
          Drop[globaltimes, -1];
globalwalk2 = Drop[globalwalk, -1];
sigmavalues = 
        Map[(1/( sigma2star[#])) &, 
          globalwalk2];
valuesA = 
        FoldList[Plus, 0, 
          sigmavalues*timesdifferences];
globalwalk3 = (Join[
              globalwalk2, {If [
                  Abs[aA - Last[globalwalk]] < Abs[bB - Last[globalwalk]], 
                  aA, bB]}] // Flatten);
\(\(ListPlot[{globaltimes, valuesA} // Transpose, 
        PlotJoined -> True, 
        PlotRange -> 
          All, 
TextStyle -> {FontFamily -> "Tahoma", 
            FontSize -> 
              14}, 
PlotLabel -> \*""Process (Subscript[A, t])""];\)\)

3.1.3.3 The population size (Subscript[N, t]), 0 ≤ t ≤ 90

population = 
      Interpolation[{valuesA, globalwalk3} // Transpose]

\(\(Plot[
        population[t], {t, 0, 
          90}, 
TextStyle -> {FontFamily -> "Tahoma", 
            FontSize -> 
              14}, 
PlotLabel -> "Size of population"];\)\)

3.2 The stochastic model with vaccination

3. 2. 1 Numerical example

We use the Kendall's parameter β[z] introduced in 2.2.1.

\(β[
        z] = \(2  β\_0\)\/\(\((1 - \(z\ β\_0\)\/γ)\) + \((1 - \(z\ β\_0\)\/γ)\)\^\(-1\)\), \ 
    0 ≤ z < \(\(γ\/β\_0\)\(.\)\)\)

In the stochastic case, the system of four differential equations with random parameters (depending on the simulated population size process (Subscript[N, t])) must be numerically solved.

We choose (Subscript[β, 0]) = 0.6, γ = 0.2, (Subscript[N, 0])=1, (Subscript[w, 0]) = 0.1 (Subscript[N, 0]), (Subscript[x, 0]) =  0.999 ((Subscript[N, 0]) - (Subscript[w, 0])), (Subscript[y, 0]) =  0.001 ((Subscript[N, 0]) - (Subscript[w, 0])), (Subscript[z, 0]) =  0,v  = 0.002.

Clear[x, y, z, w];
n0 = 1; 
    x0init = 0.999; β0 = 0.6; γ0 = 0.2; t∞ = 150; 
    v = 0.002; 
    pv = 0.1;
solsys2 = 
        NDSolve[{x'[
                t] == -population[t]/population[
                          0] β0 (2/((1 - ((z[t] β0)/γ)0) + (((1 - ((z[t] β0)/γ)0))^(-1)))) x[t] y[t] - 
                v x[t], y'[
                t] == population[t]/population[
                        0] β0 (2/((1 - ((z[t] β0)/γ)0) + (((1 - ((z[t] β0)/γ)0))^(-1)))) 
                  x[t] y[t] - γ0 y[t], z'[
                t] == γ0 y[t], w'[t] == v x[t], 
            x[0] == x0init*(1 - pv), w[0] == x0init*pv, 
            y[0] == 1 - x[0] - w[0], z[0] == 0}, {x, y, z, 
            w}, {t, 0, t∞}];
x = (x /. 
            solsys2)[[1]]; 
    y = (y /. solsys2)[[1]]; 
    z = (z /. solsys2)[[1]]; 
    w = (w /. 
            solsys2)[[1]];

The plots below show:

the number of susceptibles
the number of infectives
the number of removals
the number of vaccinated
in time interval 35 — 60 days since the infection started.

p1 = 
        Plot[{x[t]*population[t]} // Evaluate, {t, 35, 60}, 
          PlotRange ->  
            All  ];
p2 = 
        Plot[{y[t]*population[t]} // Evaluate, {t, 35, 60}, 
          PlotRange ->  All, PlotStyle -> {RGBColor[1, 0, 0]}, 
          DisplayFunction -> Identity];
p3 = 
        Plot[{z[t]*population[t]} // Evaluate, {t, 35, 60}, 
          PlotRange ->  All, PlotStyle -> {RGBColor[0, 1, 0]}, 
          DisplayFunction -> Identity];
p4 = 
        Plot[{w[t]*population[t]} // Evaluate, {t, 35, 60}, 
          PlotRange ->  All, PlotStyle -> {RGBColor[0, 0, 1]}, 
          DisplayFunction -> Identity];

Susceptibles:

Show[p1, p2, p3, p4, 
        DisplayFunction -> $DisplayFunction];

\(\(Show[p2, p4, 
        DisplayFunction -> $DisplayFunction, 
TextStyle -> {FontFamily -> "Tahoma", 
            FontSize -> 14}, 
        PlotLabel -> "Vaccinated and Infectives"];\)\)

\(\(Show[p4, 
        DisplayFunction -> $DisplayFunction, 
TextStyle -> {FontFamily -> "Tahoma", 
            FontSize -> 14}, 
        PlotLabel -> "Vaccinated "];\)\)

4 Acknowledgment

The research has been supported by Grant of Czech Ministry of Education MSM 0021620839.

5 References

[1] Abell, M. L. and Braselton, J. P.: Differential Equations with Mathematica. 3rd edition Elsevier Academic Press, Amsterdam, 2004.
[2] Daley, D. J. and Gani, J.: Epidemic Modelling; An Introduction. Cambridge University Press, Cambridge, 1999.
[3] Hlubinka, J. and Hurt, J.: A stochastic version of Kermack and Mc-Kendrick model of epidemics. In: ROBUST 2006 (Antoch, J. and Dohnal, G., eds.), JČMF, Praha, 2006. 99 - 108.
[4] Kermack, W. O. and McKendrick, A. G.: A contribution to the mathematical theory of epidemics. Proc. Roy. Soc. London A 115, 700 - 721.
[5] Kendall, D. G.: Deterministic and stochastic epidemics in closed population. In: Proc. Third Berkeley Symp. Math. Statist. Probab. 4, Univ. of California Press, Berkeley. 149 - 165.