SIR disease in Heterogeneous Populations

4.7. SIR disease in Heterogeneous Populations#

We now consider an SIR infectious disease spreading in a large population made up of heterogeneous individuals having varying contact rates. That is, some individuals have more interactions per unit time than others.

The entire population consists of $N$ individuals. Each individual has a contact rate $k$, which we will assume is an integer. Additionally we assume that the only changes that happen during the outbreak is infection and recovery. Individuals with a given contact rate maintain that contact rate for the duration of the outbreak.

4.7.1. Mathematical setup#

We will divide the susceptible, infected, and recovered individuals by their contact rate:

\[\begin{align*} S &= S_0 + S_1 + \cdots\\ I &= I_0 + I_1 + \cdots\\ R &= R_0 + R_1 + \cdots \end{align*}\]

Because individuals never change contact rate, the total number of people with each contact rate can be written as

\[ N_k = S_k+I_k + R_k \]

and the total population $N$ satisfies

\[\begin{align*} N &= N_0 + N_1 + \cdots\\ &= S+I+R \end{align*}\]

We take $p_k$ to be the proportion of the population with a given $k$:

\[ p_k = N_k/N \]

We introduce the probabiilty generating function

\[ \psi(x) = \sum_k p_k x^k \]

Then

\[\begin{align*} \sum_{k} k N_k &= \sum_k k p_k N\\ &= N \sum_k kp_k\\ &= N \sum_k k p_k 1^{k-1}\\ &= N \psi'(1) \end{align*}\]

where $\psi'(1)= \sum_k kp_k$ is the average contact rate.

Each individual in $I_{\hat{k}}$ transmits at rate $\beta \hat{k}$, which is proportional to its contact rate. So using superposition, the total rate of transmissions is

\[ \text{combined transmission rate} = \beta \sum_{\hat{k}}\hat{k} I_{\hat{k}} \]

When a transmission occurs, the recipient is chosen randomly from the population with probability proportional to its contact rate. So the probability of choosing a specific individual with contact rate $k$ is $\frac{k}{\sum_{\hat{k}} \hat{k}N_{\hat{k}}}= k/N\psi'(1)$. The probability that the recipient is both susceptible and has contact rate $k$ is the product of the overall transmission rate with the probability the recipient is one of the $S_k$ individuals with contact rate $k$:

\[ \text{combined rate of infection of individuals in $S_k$} = \beta\sum_{\hat{k}}\hat{k} I_{\hat{k}} \frac{k S_k}{N\psi'(1)} \]

For simplicity of notation we now introduce

\[ \pi_I = \frac{ \sum_{\hat{k}} \hat{k} I_k}{N\psi'(1)} \]

which is the probability that a partner chosen randomly from the population is infected. So

\[ \text{combined rate of infection of individuals in $S_k$} = \beta \pi_I k S_k \]

The per-individual recovery rate is $\gamma$, so the combined rate at which infected individuals with contact rate $k$ recover is $\gamma I_k$.

We now wave hands slightly and say that

if the numbers are large enough, it is reasonable to assume that noise due to stochastic effects is small compared to the expected change of the variables.

Then we can neglect stochastic effects and conclude that

\[\begin{align*} \frac{d}{dt} S_k &= - \beta \pi_I k S_k\\ \frac{d}{dt} I_k &= \beta \pi_I k S_k - \gamma I_k\\ \frac{d}{dt} R_k &= \gamma I_k\\ \pi_I &= \frac{\sum_k k I_k}{N\psi'(1)} \end{align*}\]

4.7.2. Reducing the number of equations#

4.7.2.1. Naive approach#

The current system of equations is infinite. One way to attempt to solve them is to truncate the system, ignoring values of $k$ larger than some threshold $K$.

However, looking at the original biological problem, the individuals with high degree get infected sooner and transmit more. So at early times, the probability of being infected is proportional to $k$ and the contribution to new infections is also proportional to $k$. Thus these individuals play a role in the dynamics that is proportional to $k^2$. Disregarding the highest degree individuals is only reasonable if $\sum_{k>K} k^2 p_k$ is very small. So typically we must choose $K$ to be quite large and in some cases, (e.g., where $p_k$ decays like $1/k^2$) the impact of large values of $k$ is never negligible.

4.7.2.2. Improved approach#

We can avoid truncating the system at all, and still arrive at a small number of equations. We start with the equations for $dS_k/dt$ and rewrite them

\[ \frac{d}{dt} S_k + \beta k \pi_I S_k = 0 \]

We multiply this equation by $Q(t)= e^{\beta k\int_0^t \pi_I(\tau)\, d\tau }$. If you have seen integrating factors before, this is how $Q$ this factor was chosen. If you have not seen integrating factors before, the key reason for this choice is that we are chosing $Q(t)$ so that it’s derivative will satisfy $dQ/dt = \beta k \pi_I Q(t)$. Thus

\[\begin{align*} Q(t) \frac{d}{dt} S_k(t) + Q(t) \beta k \pi_I S_k &=0\\ Q(t) \frac{d}{dt} S_k(t) + \left(\frac{d}{dt} Q(t)\right) S_k(t) &= 0\\ \frac{d}{dt} \left(Q(t) S_k(t)\right) &=0 \end{align*}\]

using the product rule in reverse in the final line. The choice of $Q$ made the last step possible. Integrating both sides

\[\begin{align*} Q(t) S_k(t) &= C\\ S_k(t) &= C\frac{1}{Q(t)}\\ &= C e^{-k \beta \int_0^t \pi_I(\tau) \, d\tau}\\ &= C \theta(t)^k \end{align*}\]

where we set $\theta(t) = e^{-\beta \int_0^t \pi_I(\tau) \, d\tau}$ and note that $\theta(0)=1$. Then

\[ S_k(t) = S_k(0) \theta(t)^k \]

If we assume that $S_k(0) = (1-\rho) N_k$ which corresponds to infecting a proportion $\rho$ chosen uniformly at random from the population at time $0$, we conclude

\[ S_k(t) = (1-\rho) N_k \theta(t)^k = N (1-\rho) p_k \theta(t)^k \]

So

\[\begin{align*} S(t) &= \sum_k S_k(t)\\ &= (1-\rho) N\sum_k p_k \theta(t)^k\\ &= (1-\rho) N \psi(\theta(t)) \end{align*}\]

We also find

\[\begin{align*} \frac{d}{dt} \theta(t)&= - \beta \pi_I(t)\theta(t)\\ \frac{d}{dt} I(t) &= -\frac{d}{dt} S(t) - \gamma I(t)\\ &= -(1-\rho) N \psi'(\theta(t)) \theta'(t) - \gamma I(t)\\ &= -(1-\rho) N \psi'(\theta(t)) \left(-\beta \pi_I(t)\right) \theta(t) - \gamma I(t)\\ &= \beta(1-\rho)N\pi_I(t)\theta(t)\psi'(\theta(t))- \gamma I(t)\\ \frac{d}{dt} R(t) &= \gamma I(t)\\ \pi_I(t) &= \frac{\sum_k k I_k}{N\psi'(1)} \end{align*}\]

There are a few options that will save us from trying to calculate the infinite sum $\sum_k k I_k$ at each time. One of these is to write down a differential equation for $\pi_I(t)$:

\[\begin{align*} \frac{d}{dt} \pi_I(t) &= \frac{\sum_k k \frac{d I_k}{dt}}{N\psi'(1)}\\ &= \frac{\sum_k k (-\frac{d S_k}{dt} - \gamma I_k)}{N\psi'(1)}\\ &= \frac{-\sum_k k N (1-\rho) p_k k \theta(t)^{k-1}\theta'(t)}{N\psi'(1)} - \gamma \frac{\sum_k k I_k}{N\psi'(1)}\\ &= \frac{\beta \pi_I N(1-\rho)\sum_k k^2p_k \theta^k}{N\psi'(1)} - \gamma \pi_I\\ &= \beta \pi_I (1-\rho)\frac{ \theta \frac{d}{d\theta} (\theta\psi'(\theta))}{\psi'(1)} - \gamma \pi_I \end{align*}\]

So our final system is

\[\begin{align*} \frac{d}{dt} I &= \beta(1-\rho)N\pi_I\theta\psi'(\theta)- \gamma I\\ \frac{d}{dt} \pi_I &= \beta \pi_I (1-\rho)\frac{ \theta \frac{d}{d\theta} (\theta\psi'(\theta))}{\psi'(1)} - \gamma \pi_I\\ \frac{d}{dt} \theta &= - \beta \pi_I\theta\\ S &= (1-\rho)N\psi(\theta)\\ R &= N-I-S \end{align*}\]

Where we’ve used $N=S+I+R$ to remove the $dR/dt$ equation. The initial conditions are $\theta(0)=1$, $I(0) = \rho N$ and $\pi_I(0) = \rho$ which comes from choosing a fraction $\rho$ uniformly at random to infect at time $0$.