Simple epidemic models

-

Epidemiology is the study of the spread of diseases in a population. Many students seem interested in this kind of problem due to the recent Covid-19 pandemic. Here, we develop a few simple models of the epidemic.

First, some terminology.

  • Susceptibles: Those who might succumb to the disease.
  • Infectives: Those with the disease can spread it among the susceptibles.
  • Immunes: Those who are immune to the disease. This category may include the dead and isolated.
  • Latent period: A period between the infection and onset of the disease.
  • Infectious period: A period during which the individual remains infective.
In this post, we only deal with simple epidemic models in which there are only susceptibles and infectives, and individuals do not die.

Continuous-time epidemic model (with no recovery)

Let us derive a simple continuous-time epidemic model based on the following assumptions.

  1. The population size is fixed at \(n_0 + 1\). No birth, no death.
  2. The population contains only susceptibles and infectives.
  3. At the time \(t=0\), there are \(n_0\) susceptibles and one infective. (The first introduction of the disease to the population.)
  4. \(S(t)\) is the random variable representing the number of susceptibles, not infected, at time \(t\). Thus, there are \(n_0 + 1 - S(t)\) infectives at time \(t\).
  5. The likelihood of an infection occurring in a short time interval \(\delta t\) is \[\beta S(t)[n_0 + 1 - S(t)]\delta t\] where \(\beta\) is the contact rate.
This process may be regarded as a chemical reaction between two chemical species: \(S + I \to I\) with reaction rate \(\beta\). Note that we do not consider the recovery of the infected (\(I \to S\)) so that the disease can only spread in this model.

Let \(p_n(t) = \Pr(S(t) = n), n=0, 1,\cdots, n_0\) be the probability that there are \(n\) susceptibles at time \(t\). The partition theorem reads
\[p_{n}(t + \delta t) = \beta(n+1)(n_0 - n)\delta t p_{n+1}(t) + [1 - \beta n(n_0+1 -n)\delta t]p_{n}(t) + o(\delta t)\]
for \(0 \leq n \leq n_0 - 1\), and
\[p_{n_0}(t + \delta t) = (1 - \beta n_0\delta t)p_{n_0}(t) + o(\delta t)\]
for \(n = n_0\). From these, we can derive the following differential-difference equation.
\[\frac{\mathrm{d}p_{n}(t)}{\mathrm{d}t} = \beta(n+1)(n_0 - n)p_{n+1}(t) - \beta n(n_0+1 -n)p_{n}(t), ~~ (0 \leq n \leq n_0)\]
with the initial condition
\[p_n(0) = \delta_{n,n_0}.\]
Unfortunately, the analytical solution to this differential-difference equation is unknown (so far). We need to resort to numerical simulations. In simulations, we discretize the time variable. Thus, we obtain a discrete-time model similar to the one introduced in the following section.

Discrete-time epidemic model (with recovery)

Next, let us model epidemics using a discrete-time Markov chain.
The population size is fixed to \(n\). The number of the infectives at time \(t\) is \(I(t)\), and that of susceptible is \(S(t) = n - I(t)\). Let us denote the state with \(i\) infectives by \(E_i\). Assuming only one infection or one recovery at most at each step, possible transitions are
  1. \(E_i \to E_{i-1}, E_i, E_{i+1}\) for \(i= 1, 2, \cdots, n-1\),
  2. \(E_n \to E_n, E_{n-1}\),
  3. \(E_0 \to E_0\) (an absorbing state).
See the figure below.
Figure DTM. The state transition diagram of the simple discrete-time epidemic model.


The transition probabilities are defined as
\[\begin{eqnarray} p_{i,i-1} &=& i\gamma\tau,~~(i = 1, 2, \cdots, n),\\ p_{i,i} &=& 1 - \left[\beta \frac{i(n-i)}{n}\tau + i\gamma\tau\right],~~(i = 0,1,2,\cdots, n),\\ p_{i,i+1} &=& \beta\frac{i(n-i)}{n}\tau, ~~(i=0, 1, 2,\cdots, n - 1),\\ p_{i,j} &=& 0, ~~(\text{otherwise}) \end{eqnarray}\]
where \(\tau\) is the unit time, \(\beta\) is the contact rate (as before), and \(\gamma\) is the recovery rate. In particular, the probability of a new infection is given by
\[\beta \frac{i (n - i)}{n}\tau.\]
This may be interpreted as follows. We are assuming one person is infected at a time. The probability of picking one susceptible among the population of \(n\) individuals is \((n-i)/n\); the probability of infection is proportional to the number of the infected \(i\), and the time \(\tau\) exposed to them. Thus, the probability of infection is proportional to \(\frac{n-i}{n}\times i\times \tau\).

Exercise. Given the above probabilities, verify that
\[\sum_{j=0}^{n}p_{i,j} = 1\]
for \(i = 0, 1, \cdots, n\). □

Now, let \(p_i(t) = \Pr(I(t) = i)\) be the probability that \(n\) individuals are infected at time \(t\). As usual, we can use the partition theorem to derive the following difference equations:
\[\begin{eqnarray} p_0(t + \tau) &=& \gamma \tau p_1(t) + p_0(t)\\ p_i(t + \tau) &=& \beta\frac{(i-1)(n-i+1)}{n}\tau p_{i-1}(t) + \gamma (i + 1)\tau p_{i+1}(t)\tag{Eq:DE0} \\ &&+ \left(1 - \left[\beta\frac{i(n-i)}{n}\tau + \gamma i \tau\right]\right)p_{i}(t), ~~ (i > 0).\tag{Eq:DEi} \end{eqnarray}\]

Again, we can perform numerical simulations based on this model. A sample path is shown below.


Figure SP. Purple: A simulated sample path of the discrete-time model with \(tau = 0.01\), \(\beta = 1\), \(\gamma = 0.5\), \(n = 100\), and \(I(0) = 2\). Green: Solution of the deterministic model (see the following section).

See also the video explaining how to simulate this process using Google Sheets.


Deterministic epidemic model

Next, we derive a deterministic model from the above discrete-time model. By rearranging Eqs. (Eq:DE0) and (Eq:DEi), and taking the limit \(\tau \to 0\), we have the following differential-difference equations:
\[\begin{eqnarray}\frac{dp_0(t)}{dt} &=& \gamma p_1(t),\\ \frac{dp_i(t)}{dt} &=& \beta\frac{(i-1)(n-i+1)}{n}p_{i-1}(t) + \gamma(i+1)p_{i+1}(t)\\ && -\left(\beta\frac{i(n-i)}{n} + \gamma i\right)p_i(t). \tag{Eq:dDEi}\end{eqnarray}\]

Exercise. Derive these! □

Let \(I(t)\) and \(S(t)\) be the random variables representing the number of the infectives and susceptibles, respectively, and let us define
\[\begin{eqnarray}x_I(t) &=& \mathbb{E}(I(t)) = \sum_{i=1}^{n}ip_i(t),\\x_S(t) &=& \mathbb{E}(S(t)) = \sum_{i=0}^{n-1}(n-i)p_i(t).\end{eqnarray}\]
Multiply both sides of Eq. (Eq:dDEi) and sum over \(i\), and we have the following equation:
\[\frac{dx_I(t)}{dt} = \frac{\beta}{n}\mathbb{E}(S(t)I(t)) - \gamma x_I(t)\]
where
\[\mathbb{E}(S(t)I(t)) = \sum_{i=1}^{n-1}i(n-i)p_i(t).\]
At this point, let us make an approximation:
\[\mathbb{E}(S(t)I(t)) \approx \mathbb{E}(S(t))\mathbb{E}(I(t)) = x_S(t)x_I(t),\]
then the above differential equation becomes
\[\frac{dx_I(t)}{dt} = \frac{\beta}{n}x_S(t)x_I(t) - \gamma x_I(t).\]
By the above definition, we have
\[x_S(t) + x_I(t) = n.\]
Using this equality, we can eliminate the variable \(x_S\) to have the deterministic differential equation for \(x_I(t)\):
\[\frac{{d}x_I}{{d}t} = x_I[(\beta - \gamma)n - {\beta}x_I]/n.\]
Here, we assume that 
\[(\beta - \gamma)n - \beta x_I > 0 \tag{Eq:BG}\]
always hold (this means that the infections keep increasing). With the initial condition \(x_I(0) = 2\), this can be readily solved to give
\[x_I(t) = \frac{2(\beta - \gamma)n}{2\beta + [(\beta - \gamma)n - 2\beta]e^{-(\beta-\gamma)t}}.\tag{Eq:xIsol}\]
In the limit \(t \to \infty\), 
\[x_I(t) \to \frac{(\beta - \gamma)n}{\beta}.\]
(From (Eq:BG) above, we have \(\beta - \gamma > 0\).)

Epidemics are often reported regarding the number of new cases in a certain interval (per day or week, for example). That is given by the derivative of Eq. (Eq:xIsol):
\[\frac{dx_I(t)}{dt} = \frac{2n(\beta - \gamma)^2[(\beta - \gamma)n - 2\beta]e^{-(\beta-\gamma)t}}{[2\beta + [(\beta - \gamma)n - 2\beta]e^{-(\beta-\gamma)t}]^2}.      \tag{Eq:EC}\]
The graph of this derivative is known as the epidemic curve. An example epidemic curve is shown below.
Figure EC. A sample epidemic curve, Eq. (Eq:EC). The parameters are the same as in Figure SP above.




Comments

Post a Comment

Popular posts from this blog

Open sets and closed sets in \(\mathbb{R}^n\)

-
Image
Open sets In \(\mathbb{R}\), we have the notion of an open interval such as \((a, b) = \{x \in \mathbb{R} | a < x < b\}\). We want to extend this idea to apply to \(\mathbb{R}^n\). We also introduce the notions of bounded sets and closed sets in \(\mathbb{R}^n\). Recall that the \(\varepsilon\)-neighbor of a point \(x\in\mathbb{R}^n\) is defined as \(N_{\varepsilon}(x) = \{y \in \mathbb{R}^n | d(x, y) < \varepsilon \}\) where \(d(x,y)\) is the distance between \(x\) and \(y\). Definition (Open set) A subset \(U\) of \(\mathbb{R}^n\) is said to be an  open set  if the following holds: \[\forall x \in U ~ \exists \delta > 0 ~ (N_{\delta}(x) \subset U).\tag{Eq:OpenSet}\] That is, for every point in an open set \(U\), we can always find an open ball centered at that point, that is included in \(U\). See the following figure. Perhaps, it is instructive to see what is  not  an open set. Negating (Eq:OpenSet), we have \[\exists x \in U ~ \forall \delta > 0 ~ (N_{\delta}(x) \not
Read more

Euclidean spaces

-
Image
We would like to study multivariate functions (i.e., functions of many variables), continuous multivariate functions in particular. To define continuity, we need a measure of "closeness" between points. One measure of closeness is the Euclidean distance. The set \(\mathbb{R}^n\) (with \(n \in \mathbb{N}\)) with the Euclidean distance function is called a Euclidean space. This is the space where our functions of interest live. The real line is a geometric representation of \(\mathbb{R}\), the set of all real numbers. That is, each \(a \in \mathbb{R}\) is represented as the point \(a\) on the real line. The coordinate plane , or the \(x\)-\(y\) plane , is a geometric representation of \(\mathbb{R}^2\), the set of all pairs of real numbers. Each pair of real numbers \((a, b)\) is visualized as the point \((a, b)\) in the plane. Remark . Recall that \(\mathbb{R}^2 = \mathbb{R}\times\mathbb{R} = \{(x, y) | x, y \in \mathbb{R}\}\) is the Cartesian product of \(\mathbb{R}\) with i
Read more

Applications of multiple integrals

-
Image
We can use multiple integrals to compute areas and volumes of various shapes. Area of a planar region Definition (Area) Let \(D\) be a bounded closed region in \(\mathbb{R}^2\). \(D\) is said to have an area if the multiple integral of the constant function 1 over \(D\), \(\iint_Ddxdy\), exists. Its value is denoted by \(\mu(D)\): \[\mu(D) = \iint_Ddxdy.\] Example . Let us calculate the area of the disk \(D = \{(x,y)\mid x^2 + y^2 \leq a^2\}\). Using the polar coordinates, \(x = r\cos\theta, y = r\sin\theta\), \(dxdy = rdrd\theta\), and the ranges of \(r\) and \(\theta\) are \([0, a]\) and \([0, 2\pi]\), respectively. Thus, \[\begin{eqnarray*} \mu(D) &=& \iint_Ddxdy\\ &=&\int_0^a\left(\int_0^{2\pi}rd\theta\right)dr\\ &=&2\pi\int_0^a rdr\\ &=&2\pi\left[\frac{r^2}{2}\right]_0^a = \pi a^2. \end{eqnarray*}\] □ Volume of a solid figure Definition (Volume) Let \(V\) be a solid figure in the \((x,y,z)\) space \(\mathbb{R}^3\). \(V\) is sai
Read more