Super-exponential growth of epidemics in networks with cliques

In this work, we will study the dynamics of disease propagation on complex networks containing fully connected sub-graphs or cliques. As previously noted in Refs. Mann and Dobson (2023); Bianconi and Dorogovtsev (2024); Li et al. (2021b); Karrer and Newman (2010), factor graphs provide a suitable representation for these networks. A factor graph is an undirected bipartite network formed by three sets: nodes (or individuals), factor nodes (or group nodes), and links that exclusively connect nodes with factor nodes. This concept is illustrated in Fig. 1. The left panel displays a network with cliques, while the right panel shows its associated factor graph. Alternatively, the network on the left can be interpreted as a projection of the factor graph on the right. The degree (or membership) of a node in the factor graph, $k_{I}$, corresponds to the number of links attached to it, while the degree of a factor node, $k_{C}$, is defined analogously. The number of nodes and factor nodes are denoted by $N_{I}$ and $N_{C}$, respectively. On the other hand, the degree distributions of nodes and factor nodes are represented by $P(k_{I})$ and $P(k_{C})$, respectively.

Several methods have been proposed to generate ensembles of random networks with cliques (or factor graphs) Alves et al. (2024); Nikolaev and Mneimneh (2023); Leung and Weitz (2016); Guillaume and Latapy (2006); Rizi et al. (2024). In this work, we will concentrate on the configuration model Guillaume and Latapy (2006); Rizi et al. (2024) which has been widely used to construct random and uncorrelated graphs with a prescribed degree distribution. As it was pointed out in Ref. Karrer and Newman (2010), sparse graphs generated by the configuration model exhibit a locally tree-like structure in the thermodynamic limit, and thanks to this property, many statistical metrics characterizing the structure of these networks can be calculated by using the generating-function technique. For factor graphs, the following two pairs of probability-generating functions (pgf’s) for nodes and factor nodes are commonly employed Newman et al. (2001):

	$\displaystyle G_{0}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}P(k_{I})x^{k_{I}},G_{1}(x)=\sum_{k_{I}=0}^{\infty}\frac{k_{I}P(k_{I})}{\langle k_{I}\rangle}x^{k_{I}-1},$		(1)
	$\displaystyle F_{0}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}P(k_{C})x^{k_{C}},F_{1}(x)=\sum_{k_{I}=0}^{\infty}\frac{k_{C}P(k_{C})}{\langle k_{C}\rangle}x^{k_{C}-1},$		(2)

where $x$ is a dummy variable, and:

1. 1.

   $\langle k_{I}\rangle=\sum k_{I}P(k_{I})$ and $\langle k_{C}\rangle=\sum k_{C}P(k_{C})$ are the expected degree of nodes and factor nodes, respectively.

2. 2.

   $G_{0}(x)$ is the pgf for the degree of a randomly chosen node/person in the factor graph.

3. 3.

   $G_{1}(x)$ is the generating function for the number of additional links a node has when it is reached by following a random link in the factor graph. This number is commonly referred to as the excess-degree.

4. 4.

   $F_{0}(x)$ is the pgf for the degree of a randomly chosen factor node.

5. 5.

   $F_{1}(x)$ is the generating function for the number of additional links a factor node has when it is reached by following a random link in the factor graph.

To make these concepts more concrete, consider a scientific collaboration network represented as a bipartite graph where authors/scientists (nodes) connect to their publications (factor nodes). In this example, $G_{0}(x)$ is the generating function for the probability of randomly selecting a scientist who has contributed to a number $k_{I}$ of publications. On the other hand, $F_{0}(x)$ is the generating function for the probability of randomly selecting a publication written by $k_{C}$ scientists.

Thanks to the locally tree-like character of these generated factor graphs, we can also calculate the degree distribution of their associated projected networks using the generating functions given in Eqs. (1) and (2). More specifically, Newman et al Newman et al. (2001) showed that, by using the "power property of generating functions", the probability generating function for the degree distribution of the projection can be expressed as:

$\displaystyle\mathcal{G}^{(1)}_{0}(x_{1})=G_{0}(F_{1}(x_{1}))=\sum_{n=0}^{\infty}a_{n}(x_{1})^{n},$

(3)

where $x_{1}$ is a dummy variable and $a_{n}=\frac{1}{n!}\frac{d^{n}\mathcal{G}^{(1)}_{0}(x_{1})}{d(x_{1})^{n}}|_{x_{1}=0}$ represents the probability that a randomly chosen node has $n$ first-nearest neighbors in the projected network. Returning to our scientific collaboration example, $\mathcal{G}^{(1)}_{0}(x_{1})$ generates the distribution of the total number of unique collaborators for a randomly selected scientist.

Alternatively, the pgf $\mathcal{G}^{(1)}_{0}(x_{1})$ given in Eq. (3) can also be understood from a branching process perspective. Imagine a randomly chosen node as the ancestor or root from which lineages branch off (see Fig.2). The factor nodes connected to this ancestor form the first generation and their number is captured by the outer function $G_{0}(\cdot)$ in Eq. (3). Each factor node then produces a number of "descendants" or nodes that form the second generation, and this number is encoded in the inner function $F_{1}(\cdot)$ in Eq. (3).

Another distribution that describes the local structure of projected networks, which can be easily computed using the method presented in Ref. Newman et al. (2001), is the excess-degree distribution foot01 . Newman et al Newman et al. (2001) showed that the pgf for the excess-degree distribution in these networks is given by:

$\displaystyle\mathcal{G}_{1}^{(1)}(x_{1})=G_{1}(F_{1}(x_{1})).$

(4)

By applying the same approach previously discussed, we can also derive the pgf’s for higher-order neighbors in projected networks. For instance, the pgf for the number of second-nearest neighbors of a randomly chosen individual is:

	$\displaystyle\mathcal{G}^{(2)}_{0}(x_{2})$	$\displaystyle=$	$\displaystyle G_{0}(F_{1}(G_{1}(F_{1}(x_{2})))),$		(5)
		$\displaystyle=$	$\displaystyle\mathcal{G}_{0}^{(1)}(\mathcal{G}_{1}^{(1)}(x_{2}))=\sum_{n=0}^{\infty}b_{n}(x_{2})^{n},$

where $x_{2}$ serves as a placeholder, and $b_{n}=\frac{1}{n!}\frac{d^{n}\mathcal{G}^{(2)}_{0}(x_{2})}{d(x_{2})^{n}}|_{x_{1}=0}$ is the probability that the chosen node has $n$ second-nearest neighbors. In Fig. 3, we show a graphical representation of $\mathcal{G}^{(2)}_{0}(x_{2})$. This pgf can also be understood in terms of a branching process (as shown in Fig.2): the outer function $\mathcal{G}_{0}^{(1)}(\cdot)$ in Eq. (5) corresponds to the first-neighbors of a randomly chosen node, and the inner function $\mathcal{G}_{1}^{(1)}(\cdot)$ corresponds to the first-neighbors of those first-neighbors.

Similarly, if we randomly choose an individual through a clique, the pgf for the number of outgoing second-neighbors is:

	$\displaystyle\mathcal{G}^{(2)}_{1}(x_{2})$	$\displaystyle=$	$\displaystyle G_{1}(F_{1}(G_{1}(F_{1}(x_{2})))),$		(6)
		$\displaystyle=$	$\displaystyle\mathcal{G}_{1}^{(1)}(\mathcal{G}_{1}^{(1)}(x_{2})).$		(7)

To give a concrete example of the pgf’s defined in this section, consider a random network where every person belongs to exactly two cliques, and every clique contains three members. For this case, we have the following generating functions: $G_{0}(x)=x^{2}$, $G_{1}(x)=x$, $F_{0}(x)=x^{3}$, $F_{1}(x)=x^{2}$. Then, by combining these functions, we get:

• •

  $\mathcal{G}^{(1)}_{0}(x_{1})=G_{0}(F_{1}(x_{1}))=(x_{1})^{4}$, which indicates that each person has four neighbors in the projected network.

• •

  $\mathcal{G}^{(2)}_{0}(x_{2})=G_{0}(F_{1}(G_{1}(F_{1}(x_{2}))))=(x_{2})^{8}$, which indicates that each person has eight second-neighbors in the projected network.

On top of the networks described earlier, we will investigate a compartmental susceptible-infected-quarantined (SIQ) model in which individuals who test positive are placed in quarantine (along with their neighbors). In our model, susceptible people are healthy but at risk of infection, infected people can transmit the disease to susceptible neighbors, and quarantined people are no longer in contact with the rest of the population and cannot transmit the disease. For this model, we redefine the parameter $f$ presented in Ref. Valdez et al. (2023) to denote the proportion of the population who has regular access to testing. More specifically, in our model, we differentiate between two types of individuals:

• •

  Untested people $(1-f)$: those without access to testing.

• •

  Tested people ($f$): those who have regular access to fast and accurate testing. If they test positive, they immediately quarantine themselves and then impose quarantine on their neighbors to prevent further spread.

For simplicity, we will assume that the people with and without access to testing are randomly distributed across the network nodes. On the other hand, we define $S_{\ell}$ and $S_{r}$ as the fractions of susceptible people with and without access to testing, respectively.

To describe the time evolution of the density of open cliques, we define $P^{t}_{\ell,r,i}$ as the fraction of open cliques at time $t$ containing:

• •

  $\ell$ susceptible members with access to testing,

• •

  $r$ susceptible individuals without access to testing, and

• •

  $i$ infected individuals without access to testing.

In Fig. 5, we show several examples of open cliques with various member configurations. For brevity, $P^{t}_{\ell,r,i}$ will be referred to as the fraction of open cliques containing $(\ell,r,i)$ members. Of course that, in the definition of $P^{t}_{\ell,r,i}$, we could also add another index to account for the number of quarantined people within an open clique (see for example in Fig. 4 case III and IV that $c_{3}$, which is an open clique, has a quarantined member at time $t+1$). However, we omit this index because quarantined individuals do not interact with other people, and therefore they are not relevant to the dynamics of disease spread.

Based on this multivariate probability distribution $P_{\ell,r,i}^{t}$, we proceed to define the following generating functions which are just generalizations of the pgf’s $F_{0}(x)$ and $F_{1}(x)$ given in Sec. II.1:

• •

  $F_{0}^{t}(x,y,z)=\sum_{\ell}\sum_{r}\sum_{i}P_{\ell,r,i}^{t}x^{\ell}y^{r}z^{i}$. This function generalizes $F_{0}(x)$ (defined in Sec. II.1), and corresponds to the generating function for the probability of randomly selecting an open clique with $(\ell,r,i)$ members. Note that $F_{0}^{t}(1,1,1)=\sum_{\ell}\sum_{r}\sum_{i}P_{\ell,r,i}^{t}\leq 1$ represents the total fraction of open cliques at time $t$.

• •

  $F_{1,\ell}^{t}(x,y,z)=\sum_{\ell}\sum_{r}\sum_{i}\frac{\ell P_{\ell,r,i}^{t}}{\kappa_{\ell}^{t}}x^{\ell-1}y^{r}z^{i}$. This function is analogous to $F_{1}(x)$ from Eq. (2) and represents the generating function for the probability of selecting an open clique via a member who has regular access to testing services. Here, $\kappa_{\ell}^{t}=\sum_{\ell}\sum_{r}\sum_{i}\ell P_{\ell,r,i}^{t}$.

• •

  $F_{1,r}^{t}(x,y,z)=\sum_{\ell}\sum_{r}\sum_{i}\frac{rP_{\ell,r,i}^{t}}{\kappa_{r}^{t}}x^{\ell}y^{r-1}z^{i}$. This function is analogous to the previous one and represents the generating function for selecting an open clique via a susceptible member without access to testing. Here, $\kappa_{r}^{t}=\sum_{\ell}\sum_{r}\sum_{i}rP_{\ell,r,i}^{t}$.

Moving forward, we now introduce the generating functions for the subpopulation of susceptible individuals who have regular access to testing. These functions, denoted $G_{0,\ell}^{t}(x)$ and $G_{1,\ell}^{t}(x)$, are analogous to the functions $G_{0}(x)$ and $G_{1}(x)$ given in Sec. II.1:

	$\displaystyle G_{0,\ell}^{t}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}S_{\ell}(k_{I},t)x^{k_{I}},$		(8)
	$\displaystyle G_{1,\ell}^{t}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}\frac{k_{I}S_{\ell}(k_{I},t)}{\langle k_{\ell}\rangle}x^{k_{I}-1},$		(9)

where $\langle k_{\ell}\rangle=\sum_{k_{I}=0}^{\infty}k_{I}S_{\ell}(k_{I},t)$, and $S_{\ell}(k_{I},t)$ represents the proportion of susceptible individuals at time $t$ who have regular access to testing and with a membership $k_{I}$ (meaning that an individual has $k_{I}$ cliques in the factor graph). Note that evaluating $G_{0,\ell}^{t}(x)$ at $x=1$, provides the total fraction of susceptible people with access to testing services at time $t$, i.e., $G_{0,\ell}^{t}(1)=S_{\ell}^{t}$.

Similarly, for the subpopulation of susceptible individuals without access to testing, we define the analogous generating functions:

	$\displaystyle G_{0,r}^{t}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}S_{r}(k_{I},t)x^{k_{I}},$		(10)
	$\displaystyle G_{1,r}^{t}(x)$	$\displaystyle=$	$\displaystyle\sum_{k_{I}=0}^{\infty}\frac{k_{I}S_{r}(k_{I},t)}{\langle k_{r}\rangle}x^{k_{I}-1},$		(11)

where $\langle k_{r}\rangle=\sum_{k_{I}=0}^{\infty}k_{I}S_{r}(k_{I},t)$, and evaluating $G_{0,r}^{t}(x)$ at $x=1$ provides the total fraction of the population without access to testing at time $t$ (i.e., $G_{0,r}^{t}(1)=S_{r}^{t}$).

Now, as in Sec. II.1, we can employ the power property of generating functions to obtain the pgf’s describing the local neighborhood of an individual. As an example, consider that, at time $t$, we randomly select a susceptible individual "$j$" who doesn’t have access to testing. The neighborhood composition of this person is captured by

$\displaystyle\mathcal{G}^{(1)}_{r}(x_{1},y_{1},z_{1})$

$\displaystyle=$

$\displaystyle G_{0,r}^{t}(F_{1,r}^{t}(x_{1},y_{1},z_{1}))=\sum_{\ell^{{\dagger}}=0}\sum_{r^{{\dagger}}=0}\sum_{i^{{\dagger}}=0}a_{\ell^{{\dagger}}r^{{\dagger}}i^{{\dagger}}}x_{1}^{\ell^{{\dagger}}}y_{1}^{r^{{\dagger}}}z_{1}^{i^{{\dagger}}},$

(12)

where $a_{\ell^{{\dagger}}r^{{\dagger}}i^{{\dagger}}}$ is the probability that the chosen person has in their neighborhood: 1) $\ell^{{\dagger}}$ susceptible individuals with access to testing services, 2) $r^{{\dagger}}$ susceptible people without access to these services, and 3) $i^{{\dagger}}$ infected people without access to these services. This pgf generalizes the probability-generating function presented in Eq. (3). An interesting property of this generating function that we will use extensively in our time evolution equations, is that by setting $z_{1}=0$, we are only considering those configurations where "$j$" has no infected first-nearest neighbors. Conversely, by setting $z_{1}=1$, we count all configurations in which "$j$" has any number of infected first-nearest neighbors. Similar interpretations can be made for $x_{1}=y_{1}=0$ and $x_{1}=y_{1}=1$.

On the other hand, if we randomly choose "$j$" through a clique, the corresponding pgf for the excess-degree distribution (discriminated by type) is given by:

$\displaystyle\mathcal{G}^{(1)}_{1,r}(x_{1},y_{1},z_{1})$

$\displaystyle=$

$\displaystyle G_{1,r}^{t}(F_{1,r}^{t}(x_{1},y_{1},z_{1})),$

(13)

which is a generalization of Eq. (4).

Analogously, the probability-generating functions that describe the neighborhood of individuals who have regular access to testing are:

	$\displaystyle\mathcal{G}^{(1)}_{\ell}(x_{1},y_{1},z_{1})$	$\displaystyle=$	$\displaystyle G_{0,\ell}^{t}(F_{1,\ell}^{t}(x_{1},y_{1},z_{1})),$		(14)
	$\displaystyle\mathcal{G}^{(1)}_{1,\ell}(x_{1},y_{1},z_{1})$	$\displaystyle=$	$\displaystyle G_{1,\ell}^{t}(F_{1,\ell}^{t}(x_{1},y_{1},z_{1})).$		(15)

Using the generating functions defined in the previous sections, we can now estimate the basic reproduction number $R_{0}$ at time $t=0$ for a scenario where the initial fraction of the infected population is microscopic. Because our SIQ model for $\beta=1$ is a special case of the SIRQ model explored in Valdez et al. (2023), the derivation of $R_{0}$ presented in that work remains valid here. Following Refs. Valdez et al. (2023); Valdez (2024), $R_{0}$ is calculated as the ratio between the number of infected second-nearest neighbors, and the number of infected first-nearest neighbors. The expression of $R_{0}$ is given by:

$\displaystyle R_{0}=\left(\frac{dF_{1}((1-f)x)}{dx}\Bigr{|}_{\begin{subarray}{c}x=1\end{subarray}}\frac{d\mathcal{G}^{(1)}_{1}(x)}{dx}\Bigr{|}_{\begin{subarray}{c}x=1\end{subarray}}\right)/\left(\frac{dF_{1}(x)}{dx}\Bigr{|}_{\begin{subarray}{c}x=1\end{subarray}}\right).$

(16)

To understand this expression, consider that we randomly choose an individual (referred to as the index-case) through a clique and assume that this person is initially infected. The denominator in Eq. (16) represents the average number of susceptible first-nearest neighbors who can be infected by this index-case. Once infected, those individuals without access to testing (represented by the term $\frac{dF_{1}((1-f)x)}{dx}\Bigr{|}_{\begin{subarray}{c}x=1\end{subarray}}$) will transmit the disease further to their own neighbors. Each of these untested individuals will, on average, transmit the disease to $\frac{d\mathcal{G}^{(1)}_{1}(x)}{dx}\Big{|}_{x=1}$ additional people. So the numerator represents the total average number of second-nearest neighbors of the index-case who will contract the disease.

Here, we will formulate the time-evolution equations for $P_{\ell,r,i}^{t}$, as well as those for $S_{r}(t,k_{I})$ and $S_{\ell}(t,k_{I})$. To this end, we will begin by listing the transitions that affect the number and composition of open cliques.

On one hand, it is important to note that any open clique containing at least one infected member at time $t$ will become a closed clique at $t+1$. Figure 4 illustrates this transition: for example, in case I, cliques $c_{1}$ and $c_{2}$ transition from open to closed because they contain an infected member who spreads the disease. Mathematically, open cliques containing $(\ell^{*},r^{*},i^{*})$ members with $i^{*}\geq 1$ will exit the compartment of open cliques during the time step $t\to t+1$, and become closed cliques.

On the other hand, for open cliques containing $(\ell^{*},r^{*},i^{*}=0)$ members at time $t$, the following transitions preserve the open status of these cliques:

• •

  Among susceptible members who have access to testing, a portion $\ell\leq\ell^{*}$ remains susceptible, while the rest $\ell^{*}-\ell$ transition to a quarantine state (because they receive quarantine orders from other cliques). The probability of this event is denoted as $p(\ell|\ell^{*})$. Figure 4-case IV gives an example of this transition, where a member of clique $c_{3}$ who has regular access to testing, is placed in quarantine during the time step $t\to t+1$.

• •

  Among susceptible members without access to testing: 1) a number $i$ contract the disease from other cliques (see for example, clique $c_{1}$ in Fig. 4-Case I), 2) a number $r$ remain susceptible, and 3) a number $r^{*}-i-r$ are quarantined during the time step $t\to t+1$. The probability of this event is denoted as $p(i,r|r^{*})$.

With these transition probabilities, we can now write the Markovian equations governing the fraction of open cliques containing $(\ell,i,s)$ members as,

$\displaystyle P_{\ell,r,i}^{t+1}=\sum_{\ell^{*}\geq\ell}\sum_{r^{*}\geq r+i}P_{\ell^{*},r^{*},0}^{t}\times p(\ell|\ell^{*})p(i,r|r^{*}),$

(17)

where $p(\ell|\ell^{*})p(i,r|r^{*})\equiv p(\ell,r,i|\ell^{*},r^{*},0)$ can be considered as a transition probability tensor. The exact analytical expressions for $p(\ell|\ell^{*})$ and $p(i,r|r^{*})$ are given by,

	$\displaystyle p(\ell|\ell^{*})$	$\displaystyle=$	$\displaystyle\binom{\ell^{*}}{\ell}\Phi^{\ell}\times(\mathcal{G}^{(1)}_{\ell,1}(1,0,1)-\Phi)^{\ell^{*}-\ell},$		(18)
	$\displaystyle p(i,r|r^{*})$	$\displaystyle=$	$\displaystyle\binom{r^{*}}{i,r,r^{*}-i-r}\sigma^{i}\Psi^{r}\times(1-\sigma-\Psi)^{r^{*}-i-r},$		(19)

with,

	$\displaystyle\Phi$	$\displaystyle=$	$\displaystyle\mathcal{G}^{(1)}_{1,\ell}[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0],$		(20)
	$\displaystyle\Psi$	$\displaystyle=$	$\displaystyle\mathcal{G}^{(1)}_{1,r}[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0],$		(21)
	$\displaystyle\sigma$	$\displaystyle=$	$\displaystyle G_{1,r}^{t}[F_{1,r}^{t}(0,1,1)-F_{1,r}^{t}(0,1,0)+F_{1,r}^{t}(\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0)]-\Psi,$		(22)

and their derivation is explained in detail in the Supplemental Material. An important observation is that unlike random networks without cliques, where the dynamics can be described by a single time-evolution equation Miller (2011), in our work, instead, the dynamics is governed by the system of equations given in (17) which contains an order of $\mathcal{O}((k_{C,max})^{3})$ equations, where $k_{C,max}$ is the size of the largest clique. While this clearly increases the complexity of the model, our system of equations will allow us to accurately describe the behavior of our SIQ model on networks with cliques in the thermodynamic limit.

Having computed $P_{\ell,r,i}^{t+1}$, we then proceed to calculate the new proportions of susceptible individuals with membership $k_{I}$. Specifically, for susceptible individuals without access to testing services, the value of $S_{r}(t+1,k_{I})$ is determined as follows:

$\displaystyle S_{r}(t+1,k_{I})$

$\displaystyle=$

$\displaystyle S_{r}(t,k_{I})\left(F_{1,r}^{t}\left[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0\right]\right)^{k_{I}},$

(23)

where $F_{1,r}^{t}\left[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0\right]$ is the probability that a clique remains open at time $t+1$, as explained in more detail in the Supplemental Material. Similarly, for a susceptible individual with access to testing services, $S_{\ell}(t+1,k_{I})$ is calculated as,

$\displaystyle S_{\ell}(t+1,k_{I})$

$\displaystyle=$

$\displaystyle S_{\ell}(t,k_{I})\left(F_{1,\ell}^{t}\left[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0\right]\right)^{k_{I}}.$

(24)

Once $P_{\ell,i,s}^{t+1}$, $S_{\ell}(t+1,k_{I})$ and $S_{r}(t+1,k_{I})$ are calculated, we proceed to compute other epidemiological quantities that describe the disease spread in our SIQ model. Particularly, we will focus on: 1) the total fraction of susceptible people, denoted as $S(t+1)$ and 2) the total fraction of people who contracted the disease during the time step $t\to t+1$, denoted by $I_{new}(t+1)$.

On one hand, $S(t+1)$ is obtained by simply summing all susceptible individuals with and without access to testing:

$\displaystyle S(t+1)$

$\displaystyle=$

$\displaystyle\sum_{k_{I}}S_{\ell}(t+1,k_{I})+S_{r}(t+1,k_{I}).$

(25)

On the other hand, $I_{new}(t+1)$ is calculated using the following equation:

	$\displaystyle I_{new}(t+1)$	$\displaystyle=$	$\displaystyle\mathcal{G}^{(1)}_{\ell}\left[1,1,1\right]-\mathcal{G}^{(1)}_{\ell}\left[1,1,0\right]+$		(26)
			$\displaystyle+\mathcal{G}^{(1)}_{r}\left[1,1,1\right]-\mathcal{G}^{(1)}_{r}\left[1,1,0\right],$

where

• •

  $\mathcal{G}^{(1)}_{\ell}\left[1,1,1\right]-\mathcal{G}^{(1)}_{\ell}\left[1,1,0\right]$ represents the fraction of regularly tested susceptible individuals at time $t$ with at least one infected first-neighbor,

• •

  $\mathcal{G}^{(1)}_{r}\left[1,1,1\right]-\mathcal{G}^{(1)}_{r}\left[1,1,0\right]$ represents the fraction of susceptible people without access to testing services and with at least one infected first-neighbor.

After computing $S(t+1)$ and $I_{new}(t+1)$, we proceed to update all the generating functions defined in the previous section that depend on $P_{\ell,i,s}^{t}$, $S_{\ell}(t,k_{I})$ and $S_{r}(t,k_{I})$, and finally, time is increased by 1.

In this section, we will numerically compare our theoretical predictions with the results of numerical simulations under the same initial conditions and parameter settings. In particular, we will focus here on networks with cliques where the degree distributions $P(k_{C})$ and $P(k_{I})$ follow a truncated Poisson function,

$$Pois(\lambda,k_{min},k_{max})=\begin{cases}c\frac{\lambda^{k}\exp(-\lambda)}{k!},&\text{if }k_{min}\leq k\leq k_{max}\\ 0,&\text{otherwise}\end{cases}$$

where $c$ is a normalization constant. Specifically, we will use $P(k_{C})\sim Pois(7,2,20)$ and $P(k_{I})\sim Pois(3,1,20)$ (hereafter referred to as "ER₇-ER₃" networks for brevity). On the other hand, as initial conditions, we assume that a fraction $\epsilon=10^{-6}$ of the population without access to testing is infected, while the rest of the population is susceptible.

In Fig. 6, we display $I_{tot}$ at the final stage as a function of $f$. Here, $I_{tot}$ represents the fraction of the population that has ever been infected and mathematically is computed by summing the number of new infections at each time step: $I_{tot}=\sum_{t=0}^{\infty}I_{new}(t)$. From this figure, we can see that the agreement between theory and simulations is excellent. As intuitively expected, both theory and simulations show that as more people have access to testing (higher $f$), the overall fraction of the infected population ($I_{tot}$) decreases, or in other words, it becomes easier to contain the spread of the disease.

In addition, from the figure we observe that $I_{tot}$ is rapidly suppressed around $f_{c}\approx 0.369$. This sharp behavior is consistent with earlier observations made in our previous study of the SIRQ model Valdez et al. (2023), where we noticed from stochastic simulations an abrupt transition around $f=f_{c}$. However, it is important to note that in the SIRQ model, we did not investigate the behavior of epidemics near this threshold due to the high computational cost of stochastic simulations.

Now that we have the time evolution equations for our SIQ model, we can more precisely investigate this transition by numerically integrating our equations and observing how the system behaves as $f$ approaches $f_{c}$. In contrast to hybrid transitions which require a supercritical behavior with diverging slope of the order parameter D’Souza et al. (2019), from the inset of Fig. 6, we can see that in the limit $f\to f_{c}$, $I_{tot}$ behaves as a linear function $I_{tot}=a_{i}-b_{i}(f-f_{c})$ where $a_{i}$ and $b_{i}$ are fitting parameters. Consequently, the absence of diverging slope in $I_{tot}$ indicates that our SIQ model does not undergo a hybrid transition, but rather a "Type II" abrupt phase transition (see Ref. D’Souza et al. (2019)). We also provide additional results in the Supplemental Material to show that this observation holds true for a variety of network configurations as well.

Having established the agreement between our theoretical framework and stochastic simulations for the final stage, we now turn our attention to the temporal dynamics of the SIQ model.

Figures 7a-b display the time evolution of $I_{new}$ and $S$ for $f=0.1$ and $f=0.3$ on ER₇-ER₃ networks. From these figures, one can observe that the theoretical predictions are in excellent agreement with our numerical simulation results. Additionally, in the Supplemental Material, we show that our equations also accurately predict the time evolution of the SIQ model for other network topologies.

Moving forward, we will now explore the temporal evolution of infections across different values of $f$ in more detail. From Figs. 8a-c, we notice that $I_{new}(t)$ initially increases over time as an exponential function. Moreover, the exponential growth rate $\alpha$ obeys the relationship $\alpha=\ln(R_{0})$, where $R_{0}$ is the basic reproduction number defined in Eq. (16). This relationship can be understood by examining the early stages of the epidemic, when the vast majority of the population is still susceptible. During this period, the transmission process can be approximated by a Galton-Watson branching process Spouge (2019) where each infected individual transmits the disease, on average, to $R_{0}$ members of the population. Because the infected population is still small at this point, we can safely ignore the probability of a susceptible person receiving the infection from multiple sources simultaneously. As a result, during this period, the number of new infections grows as $(R_{0})^{t}$, which can be expressed equivalently as an exponential function: $(R_{0})^{t}=\exp(\alpha t)$ with $\alpha=\ln(R_{0})$.

After this initial regime where $I_{new}$ increases exponentially, we observe that the growth of $I_{new}(t)$ gradually slows down. As illustrated in Figs. 8a and b, for values of $f$ far from the critical threshold $f_{c}\cong 0.369$, $I_{new}(t)$ follows the typical trajectory observed in many epidemic models: there is an initial exponential growth, and then a well-defined peak and a subsequent decline. This behavior is also confirmed by the effective reproduction number $R_{t}$ (estimated as $R_{t}=I_{\text{new}}(t)/I_{\text{new}}(t-1)$). As shown in the insets of Figs. 8a-b, $R_{t}$ decreases steadily, which demonstrates that for $f$ values far from $f_{c}$, the propagation starts to slow down from the very beginning of the epidemic.

However, as $f$ approaches its critical value $f_{c}$, the epidemic dynamics deviates from this typical growth pattern. As illustrated in Fig. 8c, following an initial exponential increase (with a rate $\alpha=\ln(R_{0})$), we can see that the fraction of new infections accelerates even further. This accelerated growth is also reflected in the time evolution of $R_{t}$, as shown in the inset of Fig. 8c where we can observe that $R_{t}$ increases during the early stages of the outbreak.

One possible explanation for this phenomenon is that our SIQ model becomes less effective over time because healthy individuals with regular access to testing are gradually removed from the network. A concrete example of this process is shown in Fig. 4, case IV, where we can see that during the time step $t\to t+1$, clique $c_{3}$ loses its only member who can monitor and respond to infections. As a result of this event, subsequent infections within $c_{3}$ will go undetected, allowing the disease to spread more freely through the network. In other words, this example shows that our quarantine strategy has a self-undermining effect which gradually allows the disease to spread more easily. However, it is important to note that this mechanism alone is not enough to explain the super-exponential growth because, as we will see in Fig. 8d, networks without cliques do not exhibit this phenomenon.

In order to investigate the role of cliques in the emergence of this super-exponential growth, we will now explore the time evolution of the effective reproduction number $R_{t}$ in the $\langle k_{c}\rangle$-$f$ parameter space, where $\langle k_{c}\rangle$ is the average clique size of an "ER_λ-ER₃" network. Specifically, for each point in this parameter space (where $R_{0}>1$), we numerically integrate our time evolution equations for ten time steps and analyze the behavior of $R_{t}$ from $t=3$ to $t=10$ foo . If $R_{t}$ monotonically decreases during this period, we will say that $I_{new}(t)$ exhibits the standard exponential growth. On the other hand, if $I_{new}(t)$ does not decrease monotonically, this is an indication that the epidemic growth is faster than exponential. Our results shown in Fig. 8d reveal that this accelerated behavior becomes more apparent near the critical curve where $R_{0}=1$, and consequently, this suggests that the emergence of accelerated transmission dynamics could serve as an indication that the system is close to the transition point. On the other hand, from the figure we observe that the phenomenon of faster-than-exponential growth seems to be absent in networks with smaller clique sizes ($\langle k_{c}\rangle$). Therefore, our results indicate that larger cliques play a crucial role in facilitating the acceleration of disease spread. In the Supplemental Material, we show that these observations qualitatively hold for other networks with cliques with different distributions.

We assume that the fraction of people with and without access to testing are $f$ and $1-f$, respectively, and they are randomly distributed across the network nodes. Additionally, in order to reduce the number of equations that govern the dynamics between $t=0$ and $t=1$, we introduce the assumption that only a fraction $\epsilon$ of the people without access to testing services are initially infected while the rest of the population is susceptible. Thus, at $t=0$ we have that:

• •

  the total fraction of susceptible people with access to testing is $S_{\ell}(t=0)=f$, while for those without access, it is $S_{r}(t=0)=(1-f)(1-\epsilon)$,

• •

  the total fraction of infected people with access to testing is $I_{\ell}(t=0)=0$, while for those without access, it is $I_{r}(t=0)=\epsilon f$,

• •

  the fraction of susceptible people with membership $k_{I}$ and with access to testing is $S_{\ell}(t=0,k_{I})=P(k_{I})f$, and

• •

  the fraction of susceptible people with membership $k_{I}$ and without access to testing is $S_{r}(t=0,k_{I})=P(k_{I})(1-f)(1-\epsilon)$.

On the other hand, the probability of randomly selecting a clique with ($\ell,r,i$) members at time $t=0$ is,

$\displaystyle P_{\ell,r,i}^{t=0}=\sum_{k_{C}}P(k_{C})\binom{k_{C}}{\ell,r,i}f^{\ell}[(1-\epsilon)(1-f)]^{r}(\epsilon f)^{i}\delta_{\ell+r+i,k_{C}},$

(S1)

where $\delta$ is the Kronecker delta which ensures that the total number of members in the clique is equal to $k_{C}$.

This section explains how to calculate $p(\ell|\ell^{*})$ which is the probability of the following events occurring in an open clique with $\ell^{*}$ susceptible members with access to testing at time $t$: 1) exactly $\ell^{*}-\ell$ members will move to the $Q$ compartment during the time step $t\to t+1$, 2) $\ell$ members will remain susceptible. As discussed in the main text, Sec. III.3, these two events preserve a clique’s open status. Now, we will calculate the probabilities of these events below.

Here, we will calculate $p(i,r|r^{*})$ which is the probability of the following events occurring in an open clique with $r^{*}$ susceptible members without access to testing at time $t$: 1) exactly $i$ members will contract the disease from other cliques during the time step $t\to t+1$, 2) $r$ members will remain susceptible, and 3) $r^{*}-i-r$ will move to the $Q$ compartment.

Consider a randomly chosen clique denoted by “$C$”. For each susceptible member "$j$" without access to testing, three possible transitions exist:

1. 1.

   “$j$” becomes infected with probability $\sigma$.

2. 2.

   “$j$” remains susceptible with probability $\Psi$.

3. 3.

   “$j$” is quarantined with probability $1-\sigma-\Psi$.

Assuming that susceptible members without access to testing are independent of one another, we can express the transition probability $p(i,r|r^{*})$ as:

$\displaystyle p(i,r|r^{*})$

$\displaystyle=$

$\displaystyle\binom{r^{*}}{i,r,r^{*}-i-r}\sigma^{i}\Psi^{r}\times(1-\sigma-\Psi)^{r^{*}-i-r}.$

(S5)

In this equation, $\Psi$ and $\sigma$ are given by:

• •

  $\Psi=\mathcal{G}^{(1)}_{1,r}[\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0]$, which has a similar interpretation to Eq. (S3).

• •

  $\sigma=G_{1,r}^{t}[F_{1,r}^{t}(0,1,1)-F_{1,r}^{t}(0,1,0)+F_{1,r}^{t}(\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0)]-\Psi$, which is obtained through a detailed enumeration of configurations where "$j$" becomes infected during the time step $t\to t+1$. The first term represents the probability that "$j$" either remains susceptible or becomes infected while the second term, $\Psi$, is the probability that "$j$" remains susceptible. Now, focusing on the first term, the expression:

  • –

    $F_{1,r}^{t}(0,1,1)-F_{1,r}^{t}(0,1,0)$ represents the probability that an open clique (in which "$j$" is a member) has at least one infected member, and no susceptible member with access to testing.

  • –

    $F_{1,r}^{t}(\mathcal{G}^{(1)}_{1,\ell}(1,1,0),1,0)$ is the probability that an open clique (in which "$j$" is a member) has members with access to testing who remain disease-free and, therefore, do not initiate a quarantine order.