Age-Stratified COVID-19 Spread Analysis and Vaccination: A Multitype Random Network Approach

Let us consider a multitype network which consists of $M$ types of nodes, each corresponding to an age group in a population. We use $w_{i}$ to represent the fraction of the nodes of type $i\in[1,M]$. The contact from a type-$i$ node to others follows degree distribution $p_{i}(k_{1},k_{2},...,k_{M})\triangleq p_{i}(\boldsymbol{k})$, describing the joint probability for type-$i$ node to be connected with $k_{1}$ type-$1$ node, $k_{2}$ type-$2$ node, …, and $k_{M}$ type-$M$ node, where ${\boldsymbol{k}}=(k_{1},k_{2},\ldots,k_{M})$. The considered network can be generated by the following procedure: 1) generate stubs for every node following degree distribution $p_{i}(\boldsymbol{k})$, where each stub contains the information about which type of node it reaches. 2) randomly wire two matching stubs together to create an edge and repeat this process until no stubs left.

Susceptible-Infected-Removed (SIR) and Susceptible-Exposed-Infected-Removed (SEIR) compartmental models are widely used for epidemic modeling. In SEIR compartmental model, each individual can be in one of the four states, i.e., Susceptible, Exposed, Infected, or Removed. Here, to capture the salient features of COVID-19, we present a novel compartmental model, i.e., SEAHIR model, which adds two additional compartments, i.e., asymptomatic and hospitalized, to the classic SEIR model. In SEAHIR model, each individual can be in one of the six states: susceptible (S), exposed (E), symptomatic and infectious (I), asymptomatic and infectious (A), hospitalized (H), and removed (R). Both new compartments are paramount to describe the dynamics of COVID-19: the number of people in H state indicates the hospitalizations, which must be kept lower than health care capacity; patients in A state have different level of infectivity compared with symptomatic ones[27]. We assume that individuals in E state is not infectious because of low virus load, and individuals in H state are properly isolated. Individuals in I and A states are assumed to be infectious to others, where the infection rate from a type-$i$ source node to a type-$j$ node is $\lambda_{i,j}^{I}$ or $\lambda_{i,j}^{A}$ given the source node belongs to I or A state. For conciseness, we do not distinguish recovery and death in R state, but assume that an age-dependent fraction of infected people will die. Furthermore, given the fact that the cases of reinfection with COVID-19 are still extremely rare, we do not consider the transition from R state to S state. For type-$i$ nodes, the transitions among the compartments are illustrated in Fig. 2, where the symbols on the arrows denote the corresponding transition rates from one to another, which are all dependent on node type $i$ to account for the age-dependent effects.

Disease propagates from infectious nodes to its neighbors, leading to an epidemic if the epidemic size is comparable to the whole population. We employ the edge-based compartmental method to solve the equations of dynamics governing the epidemic spread over random graph[23]. The core idea of the edge-based method is to shift our attention from an individual node to the status of its neighbor reached by an edge. To study the impact of vaccination or natural immunity, we consider that a fraction of population may be already immune at the beginning of analysis: $S_{i}(\boldsymbol{k},0)$ represents the fraction of type-$i$ nodes with degree $\boldsymbol{k}$ that are initially susceptible. Besides, we use $\theta_{ji}(t)$ to denote the probability that a type-$j$ neighbor has not transmitted the disease to a type-$i$ node by time $t$ given that the type-$i$ node is susceptible at time $0$. $\theta_{ji}(t)$ can be interpreted as a state of a type-$j$ neighbor of an initially susceptible type-$i$ node. It is noted that $\theta_{ji}(0)=1$ according to the definition. Based on Fig. 2, we can construct the following equations to characterize the time-dependent epidemic process.

where $S_{i}(t)$, $A_{i}(t)$, $I_{i}(t)$, $H_{i}(t)$, $R_{i}(t)$, and $E_{i}(t)$ represent the proportions of type-$i$ nodes in the corresponding states at time $t$, respectively. From Markov chain theory, when the network size is sufficiently large, the fraction of nodes in $A_{i}$, $I_{i}$, $H_{i}$, and $R_{i}$ states can be described well by the differential equations (2)-(5) due to the flow diagram in Fig. 2. Moreover, $(\ref{d6})$ is obtained from $S_{i}(t)+A_{i}(t)+I_{i}(t)+H_{i}(t)+R_{i}(t)+E_{i}(t)=1$. For the initial conditions, we assume $A_{i}(0)=I_{i}(0)=H_{i}(0)=E_{i}(0)=0$ and $R_{i}(0)=1-S_{i}(0)$. Obviously, one can solve the above equations as long as the key probability, i.e., $\theta_{ji}(t)$, is derived. To calculate $\theta_{ji}(t)$, following the approach in [23], we break it into six parts, i.e., $\xi_{ji}^{S}(t)$, $\xi_{ji}^{E}(t)$, $\xi_{ji}^{A}(t)$, $\xi_{ji}^{I}(t)$, $\xi_{ji}^{H}(t)$, and $\xi_{ji}^{R}(t)$. Specifically, $\xi_{ji}^{S}(t)$, $\xi_{ji}^{E}(t)$, $\xi_{ji}^{A}(t)$, $\xi_{ji}^{I}(t)$, $\xi_{ji}^{H}(t)$, or $\xi_{ji}^{R}(t)$ represents the probability that the considered type-$j$ neighbor is in $S$, $E$, $A$, $I$, $H$, or $R$ state, respectively, and has not transmitted the disease to the initially susceptible type-$i$ node by time $t$, satisfying

A type-$j$ neighbor in $S$ state cannot infect the type-$i$ node, and hence $\xi_{ji}^{S}(t)$ is simply equal to the probability that the considered type-$j$ neighbor is susceptible. The degree distribution of the considered type-$j$ neighbor is given by $\frac{k_{i}p_{j}(\boldsymbol{k})}{\overline{k}_{ji}}$, where $\overline{k}_{ji}=\sum_{\boldsymbol{k}}k_{i}p_{j}(\boldsymbol{k})$ is the average degree leaving from type-$j$ node to type-$i$ node, which normalizes the probability distribution. This quantity is proportional to $k_{i}p_{j}(\boldsymbol{k})$ because type-$j$ nodes with more edges incident to type-$i$ nodes are more likely to become the neighbors of type-$i$ nodes[21]. We can obtain $\xi_{ji}^{S}(t)$ by computing the probability that the type-$j$ neighbor is initially susceptible and has not been infected by any of its neighbors, except the considered type-$i$ node, by time $t$, i.e.,

where $\delta_{il}$ is the Kronecker delta operator, with $\delta_{il}=1$ only when $i=l$, and $\delta_{il}=0$ otherwise.

As mentioned before, two states in our compartmental model, i.e., $I$ and $A$ states, are infectious. Thus, the decrease in $\theta_{ji}$ comes from two joint events: 1) the type-$j$ neighbor is $A$ or $I$ state, and 2) it transmits the disease to the type-$i$ node with infection rate $\lambda_{ji}^{A}$ or $\lambda_{ji}^{I}$, i.e.,

One can also interpret (9) in this way: there is a state $1-\theta_{ji}(t)$ (which corresponds to that the type-$j$ neighbor has transmitted the disease to the initially susceptible type-$i$ node by time $t$) receiving the flows from both $\xi_{ji}^{A}(t)$ and $\xi_{ji}^{I}(t)$.

If staying in $A$, $I$ or $H$ state, the type-$j$ neighbor transits to R state with rate $\gamma_{i}^{A}$, $\gamma_{i}^{I}$ or $\gamma_{i}^{H}$, which means

The type-$j$ neighbor progresses from I state to H state with rate $\eta_{j}$, leading to

States $\xi_{ji}^{A}(t)$ and $\xi_{ji}^{I}(t)$ receive the flows from state $\xi_{ji}^{E}(t)$. Besides, $\xi_{ji}^{A}(t)$ progresses to state $1-\theta_{ji}(t)$ and $\xi_{ji}^{R}(t)$, while $\xi_{ji}^{I}(t)$ progresses to state $1-\theta_{ji}(t)$, $\xi_{ji}^{R}(t)$, and $\xi_{ji}^{H}(t)$, yielding

In summary, one can solve $\theta_{ji}(t)$ from the following equations.

where $\xi_{ji}^{R}(0)=1-\sum_{\boldsymbol{k}}\frac{k_{i}S_{j}(\boldsymbol{k},0)p_{j}(\boldsymbol{k})}{\overline{k}_{ji}}$, $\theta_{ji}(0)=1$, and $\xi_{ji}^{A}(0)=\xi_{ji}^{I}(0)=\xi_{ji}^{H}(0)=0$. By plugging $\theta_{ji}(t)$ into (1), we can obtain the fraction of type-$i$ nodes in each compartment at given time from (1)-(6). As a result, the desired age-stratified epidemic dynamics can be obtained by solving $\mathcal{O}(M^{2})$ equations, where $M$ is the number of age groups. One can see that (1)-(6) and (14) account for considerably more population structure than a fully mixing model by only introducing marginally more complexity.

Epidemic size measures the fraction of people eventually getting infected, and epidemic probability is defined as the likelihood that the first infected patient sparks an epidemic[9]. SEAHIR model contains more compartments than traditional SIR or SEIR models, which complicates the derivations of these two key metrics. Fortunately, both metrics only depend on final state of networks. From Fig. 2, it is intuitive to see that the network progresses to an equilibrium when $t\rightarrow\infty$, where $E_{i}(\infty)=A_{i}(\infty)=I_{i}(\infty)=H_{i}(\infty)=0$. By using this fact, we use a single compartment $\mathbb{I}$ to replace all the infected states, i.e., E, A, I, H states in the original SEAHIR compartmental model. Here is our trick: although S$\mathbb{I}$R model cannot be used to capture the temporal dynamics of SEAHIR model, it can be calibrated appropriately to share the same final network state, i.e., $S_{i}(\infty)$ and $R_{i}(\infty)$, with SEAHIR model. Let $T_{ji}\in[0,1]$ denote the probability that an infected type-$j$ node ultimately transmits the disease to an initially susceptible adjacent type-$i$ node. We assume that the S$\mathbb{I}$R model is with infection rate $\hat{\lambda}_{ji}$ from type-$j$ node to type-$i$ node and transition rate $\hat{\gamma}_{ji}$ from $\mathbb{I}$ state to R state for a type-$j$ neighbor of a type-$i$ node. Since $T_{ji}$ for the S$\mathbb{I}$R model is $\frac{\hat{\lambda}_{ji}}{\hat{\lambda}_{ji}+\hat{\gamma}_{ji}}$, the S$\mathbb{I}$R model has exactly the same final network state as the SEAHIR model as long as $\frac{\hat{\lambda}_{ji}}{\hat{\lambda}_{ji}+\hat{\gamma}_{ji}}$ is set to the $T_{ji}$ in the SEAHIR model.

Let us first calculate $T_{ji}$ in the SEAHIR model, and then derive the desired metrics based on the S$\mathbb{I}$R model. An infected type-$j$ node may enter one of two infectious states, i.e., $I$ state or $A$ state, with probability $\frac{\delta_{j}}{\beta_{j}+\delta_{j}}$ or $\frac{\beta_{j}}{\beta_{j}+\delta_{j}}$. Given that the type-$j$ node in I or A state, suppose that there are two stages that it will progress to, i.e., “infecting the adjacent type-$i$ node” and “leaving current state”, with rates $\lambda_{ji}^{I}$ (or $\lambda_{ji}^{A}$) and $\eta_{j}+\gamma^{I}_{j}$ (or $\gamma^{A}_{j}$), respectively. $T_{ji}$ is the probability that the considered node enters the first stage, which equals $\frac{\lambda^{I}_{ji}}{\lambda^{I}_{ji}+(\eta_{j}+\gamma^{I}_{j})}$ if the type-$j$ node is in $I$ state, and equals $\frac{\lambda^{A}_{ji}}{\lambda^{A}_{ji}+\gamma^{A}_{j}}$ if it is in $A$ state, yielding

In S$\mathbb{I}$R model, we use $\hat{\xi}_{ji}^{S}(t)$, $\hat{\xi}_{ji}^{\mathbb{I}}(t)$, or $\hat{\xi}_{ji}^{R}(t)$ to represent the probability that a type-$j$ neighbor of an initially susceptible type-$i$ node is in $S$, $\mathbb{I}$, or $R$ state and meanwhile has not infected this type-$i$ node by time $t$, and use $\hat{\theta}_{ji}(t)$ to denote the probability that the type-$j$ neighbor has not infected the initially susceptible type-$i$ node by time $t$, satisfying $\hat{\theta}_{ji}(t)=\hat{\xi}_{ji}^{S}(t)+\hat{\xi}_{ji}^{\mathbb{I}}(t)+\hat{\xi}_{ji}^{R}(t)$. Notice that $\hat{\theta}_{ji}(t)=\theta_{ji}(t)$ only when $t=\infty$, because SEAHIR and S$\mathbb{I}$R model share the same final network state while having different temporal dynamics. By analogy with the preceding subsection, we can obtain

where $\hat{\xi}_{ji}^{R}(0)=1-\sum_{\boldsymbol{k}}\frac{k_{i}S_{j}(\boldsymbol{k},0)p_{j}(\boldsymbol{k})}{\overline{k}_{ji}}$ is the probability that the type-$j$ neighbor is initially removed (immune), and $\hat{\theta}_{ji}(0)=1$. The derivation of (16) and (18) is similar to that of (8) and (9). Analogous to the preceding subsection, probability $1-\hat{\theta}_{ji}(t)$ corresponds to that the type-$j$ neighbor has transmitted the disease to the initially susceptible type-$i$ node by time $t$. State $\hat{\xi}_{ji}^{\mathbb{I}}(t)$ transits to state $1-\hat{\theta}_{ji}(t)$ with rate $\hat{\lambda}_{ji}$ while transiting to state $\hat{\xi}_{ji}^{R}(t)$ with rate $\hat{\gamma}_{ji}$, leading to the relationship between $\hat{\xi}_{ji}^{R}(t)$ and $1-\hat{\theta}_{ji}(t)$ in (17).

Clearly, the epidemic dynamics can be governed by (18) and (19). By taking (19) into (18), we have

Now the advantage of using S$\mathbb{I}$R compartmental model becomes clearer: we arrive at (20) which is only related to $\hat{\theta}_{ji}(t)$. Since the population is closed, the epidemic will eventually go extinct, implying $\hat{\xi}_{ji}^{\mathbb{I}}(\infty)=0$. Given that $\dot{\hat{\theta}}_{ji}(\infty)=-\hat{\lambda}_{ji}\hat{\xi}_{ji}^{\mathbb{I}}(\infty)=0$ and $\theta_{ji}(\infty)=\hat{\theta}_{ji}(\infty)$, and recall that $\frac{\hat{\lambda}_{ji}}{\hat{\lambda}_{ji}+\hat{\gamma}_{ji}}=T_{ji}$, we can go back to the original SEAHIR model and obtain

where $T_{ji}$ can be calculated from (15). One can solve $\theta_{ji}(\infty)$ from the above equation. In fact, (21) can be explained in an intuitive way: a type-$i$ node has not been infected by a type-$j$ neighbor since $t=0$ either because it cannot be reached from its neighbor with probability $1-T_{ji}$, or it can be reached from its neighbor with probability $T_{ji}$ but the neighbor has not been infected since $t=0$. Thus, as $t\rightarrow\infty$, the fraction of type-$i$ nodes that have been infected since $t=0$ is given by $R_{i}(\infty)-R_{i}(0)$, i.e.,

The epidemic size is therefore expressed as

It is noted that when the population is fully susceptible, i.e., $S_{i}(\boldsymbol{k},0)=1$ for all $j$ and $\boldsymbol{k}$, $\mathcal{R}$ solved from (21) and (23) is the size of giant component in multitype networks obtained by percolation theory[21, 22]. We can compute epidemic probability in a similar way. Let $\mu_{ij}(\infty)$ denote the probability that a type-$i$ node (which is assumed to be infected) does not spark an epidemic via a type-$j$ neighbor. Analogous to (21), the following argument is true: a type-$i$ node does not ignite an epidemic via a type-$j$ node either because it cannot infect its type-$j$ neighbor, or it can infect its type-$j$ neighbor but the latter cannot spark an epidemic, which means

As a result, the probability that an infected type-$i$ node ignites an epidemic is given by

By randomly choosing a susceptible node to infect, we define the likelihood that it starts an epidemic as the epidemic probability:

where $\overline{w}_{i}=\frac{\sum_{\boldsymbol{k}}w_{i}S_{i}(\boldsymbol{k},0)}{\sum_{i}\sum_{\boldsymbol{k}}w_{i}S_{i}(\boldsymbol{k},0)}$ is the probability that the randomly chosen susceptible node is of type $i$, which reduces to $w_{i}$ in a fully susceptible population.

One of the fundamental parameters for an epidemic is its reproduction number, i.e., the average number of secondary cases caused by an infected individual. Again, since this metric is unrelated to temporal dynamics, we can derive it from S$\mathbb{I}$R model to simplify our calculation. According to the seminal work [28], the reproduction number is equivalent to the spectral radius of the next generation matrix $FV^{-1}$, where the $(m,n)$ element of matrix $F$ is the rate of new infections entering infected state $m$ caused by infected state $n$, and $(m,n)$ element of matrix $V$ is the rate at which infected state $m$ transfers to infected state $n$, assuming that the population remains near the disease-free equilibrium. In our edge-based S$\mathbb{I}$R model (18) and (19), $\hat{\xi}_{ji}^{\mathbb{I}}$ can be treated as the infected state. Therefore, we have $M^{2}$ infected states in total. Recall that there should be $\frac{\hat{\lambda}_{ji}}{\hat{\lambda}_{ji}+\hat{\gamma}_{ji}}=T_{ji}$, where $T_{ji}$ is calculated from (15). To simplify the derivations below, we further assume that $\hat{\lambda}_{ji}+\hat{\gamma}_{ji}=1$ (one can derive the same $R_{0}$ in (30) without this assumption, because $R_{0}$ is only related to $\frac{\hat{\lambda}_{ji}}{\hat{\lambda}_{ji}+\hat{\gamma}_{ji}}$). By differentiating (19) and taking (18) into it, we have

Then, to obtain the linearized subsystem for infected states about the disease-free equilibrium, we linearize (27) at the origin ($\hat{\xi}_{ji}^{\mathbb{I}}(t)=0$ and $\hat{\theta}_{ji}(t)=1$):

Then, we can construct $F$ as an $M^{2}\times M^{2}$ matrix, with $\big{(}(i-1)M+j,(j-1)M+l\big{)}$ element equal to

and construct $V$ as an $M^{2}\times M^{2}$ identity matrix. Reproduction number $R_{0}$ is therefore given by

where $\rho(\cdot)$ represents spectral radius. From Theorem 2 in [28], $R_{0}=1$ marks the epidemic threshold in the sense that the disease-free equilibrium is asymptotically stable if $R_{0}<1$, and is unstable if $R_{0}>1$. In particular, in the case of $S_{j}(\boldsymbol{k},0)=1$ for all $j$ and $\boldsymbol{k}$, $R_{0}$ becomes the basic reproduction number, i.e., the average number of secondary cases caused by an infected individual in a completely susceptible population. In this special case, the epidemic threshold $R_{0}=1$ is also in agreement with the threshold for multitype random network obtained by percolation theory[21, 22].