The rapid spread of COVID-19 has devastated the world since its inception in December 2019, leading to global economic crises and claiming hundreds of thousands of lives. As schools and businesses reopen, it is of paramount importance to asses how various safety measures may limit the spread of COVID-19. One such measure is mask-wearing, which is known to reduce the transmissibility of viruses that spread through respiratory droplets. Much of the existing work surrounding the effectiveness of mask-wearing have studied how it limits transmission between individuals [11] , [9] . However, several questions remain regarding the health of the general public. How many people must wear masks to significantly curb the spread of COVID-19? More generally, how does mask-wearing change the spreading dynamics of an epidemic?
In this paper, we provide quantitative answers to the questions above. To do so, we consider a natural generalization of the commonly-used Susceptible-Infected-Recovered (SIR) model on networks in which some individuals wear masks while others do not. We allow for different probabilities of transmission between mask-wearers and non-mask-wearers, so that an individual wearing a mask is less likely to be infected. We refer to this model as the mask model. For networks with a given degree distribution, we provide 1 yurunt@andrew.cmu.edu, 2 anirudhs@princeton.edu 3 oyagan@andrew.cmu.edu, 4 poor@princeton.edu analytical methods to accurately predict the total number of infected individuals of each type (mask-wearing and nonmask-wearing) as well as the probability that an epidemic will emerge. Technically, this is achieved by adapting techniques developed by Alexander and Day [2] as well as Eletreby, Zhuang, Carley, Yagan and Poor [8] , which were used to study a multi-strain model with mutation. Finally, we conduct extensive simulations to illustrate how mask-wearing can impact the spread of an epidemic.
Classical models of epidemics use a system of ordinary differential equations (ODEs) to describe the fraction of susceptible, infected and recovered individuals within the population (see for instance [5] ). Prior models which incorporate the effects of mask-wearing have modified the basic ODE model in various ways. Brienen et al [6] considered a simple modification in which the reproductive number of the virus, R 0 , is reduced by a multiplicative factor based on the efficacy of masks. Subsequently, Tracht et al [17] as well as Eikenberry et al [7] considered more complex generalizations of the basic ODE model in which mask-wearers and nonmask-wearers have different transmissibilities and maskwearers become non-mask-wearers at some rate, as well as vice versa. While ODE-based models are relatively simple to simulate and analyze, they are only mathematically justified under the unrealistic assumption that an infected individual can transmit the virus to any other susceptible individual in the population, regardless of location or other factors.
Our approach, on the other hand, falls under the class of network epidemic models. These models take an individuallevel view of viral spread, and studies how the structure of the contact network influences the epidemic. This provides much finer information about the epidemic, but is costly to simulate, spurring a large body of work devoted to deriving analytical predictions of epidemiological properties [15] , [14] , [12] . In particular, our work is closely related to literature on heterogenous bond percolation [3] , [10] and a multiple-strain model with mutations [2] , [8] . We elaborate on these connections in later sections.
The most basic model of network epidemics was studied by Newman [15] ; we briefly review his setup in order to provide context for the more complex models we consider in this paper. Given a prescribed degree distribution (for instance Poisson or Power law), a random contact network is generated via the configuration model [13] , [4] , [16] . Initially, a single individual (patient zero) is infected with the virus, and each neighbor of patient zero becomes infected with probability T , where T is referred to as the transmissibility of the virus. Patient zero then recovers and is no longer susceptible. The process continues as each newlyinfected vertex attempts to infect their susceptible neighbors in the same manner. The process terminates when there are no more susceptible vertices in the population.
To account for the effects of mask-wearing on viral spread, we make the following modifications to Newman's model. First, we specify m ∈ [0, 1] to be the expected fraction of individuals who wear a mask. Formally, we assign each vertex in the contact network a mask with probability m and no mask with probability 1 − m. This is done independently for each vertex. Second, we assume that the transmissibilities are heterogenous: the probability that individual u infects individual v depends on whether u and v are wearing masks. We say that a vertex is of type 1 if they wear a mask and type 2 if they do not wear a mask. We then have four parameters describing the transmissibility of the virus: T 11 , T 12 , T 21 and T 22 . The parameter T 11 is the transmissibility when u and v both wear masks, T 12 is the transmissibility when u wears a mask and v does not, etc. For brevity, we refer to this model as the Mask Model. This type of model is sometimes called heterogenous bond percolation over multi-type networks. We remark that while Allard et al [3] consider such a model in full generality, an important contribution of this paper is to study in detail the important case of mask-wearing. After the initial submission of our paper in September 2020, Lee and Zhu [10] studied the same Mask Model we propose here and derived the epidemic threshold and expected epidemic size using similar techniques as Allard et al. Here, using different techniques we additionally characterize the probability of emergence and provide extensive simulations to support our results.
In [2] , Alexander and Day proposed a multiple-strain model that accounts for mutations between strains. In their model, there are d possible strains of a virus with transmissibilities given by Q 1 , . . . , Q d . If an individual is infected with strain i, the virus may mutate into a different strain within the host. Formally, the probability that strain i mutates into strain j within a host is given by µ ij .
We next describe a mapping between the Mask Model and the multi-strain model with mutation. The key insight is that in expectation, a mask-wearing individual will have a different effective transmissibility than a non-mask-wearing individual. This will allow us to map the mask-wearing model into a two-strain model with mutation.
We begin by deriving the transmissibilities of the twostrain model. Suppose that a vertex v is infected and wears a mask. Since each neighbor wears a mask with probability m, the expected transmissibility of v is given by
Similarly, if v does not wear a mask, the transmissibility is given by
Proceeding with the analogy, the mutation probability µ 11 is the fraction of mask-wearing neighbors infected by a maskwearer. This is given by
Using the same reasoning, we can compute the other three mutation probabilities as
The advantage of this formulation is that it allows us to compute analytical predictions for the probability of emergence and epidemic size in the Mask Model using the methods of Eletreby et al [8] , which hold for the multi-strain model with mutation.
In this section, we derive analytical predictions for the probability of emergence and the expected epidemic size. One way to do so is by formulating the mask model as a multi-strain model with mutation and then leverage the analytical predictions of Eletreby et al [8] . We also compute the probability of emergence and epidemic size directly for the mask model, using methods developed by Alexander and Day [2] as well as Eletreby et al [8] .
Emergence refers to the event where the epidemic process persists over time and keeps infecting susceptible individuals. Extinction, on the other hand, is the event where the epidemic dies out in finite time. In this section, we show how to compute P 1 (resp. P 2 ) which is the probability of extinction given that patient zero wears a mask (resp. does not wear a mask). The probability of emergence can then be computed as 1 − P 1 if patient zero wears a mask and 1 − P 2 otherwise.
Our analysis follows the method of Alexander and Day [2] , who derived expressions for the probability of emergence in the multi-strain model with mutation. Suppose that a randomly chosen vertex v is patient zero and assume that v wears a mask. Let X (resp. Y ) be the number of mask-wearing (resp., non-mask-wearing) neighbors of v who are infected by v. Then, conditioned on v having k 1 susceptible mask-wearing neighbors and k 2 susceptible nonmask-wearing neighbors, X and Y are independent with X ∼ Binomial(k 1 , T 11 ) and Y ∼ Binomial(k 2 , T 12 ). Thus for s, t ∈ R,
Next, if k is the total number of susceptible neighbors, we have k 1 ∼ Binomial(k, m) and
Using the analogy with the multi-strain model with mutation (see (1)- (6)) we can equivalently write
We remark that (7) also holds in the multi-strain model with mutation [2, Section 2.2], implying the probability of emergence is identical in both models. For completeness, we describe how to compute this probability. Define γ 1 (s, t) to be the probability generating function (PGF) for the number of infections of each type (maskwearing or not) emanating from patient zero. We have
where g is the PGF of the degree distribution, i.e., g(z) = ∞ k=0 p k z k . Following the same arguments, we have γ 2 (s, t) = g(1 − Q 2 + Q 2 (sµ 21 + tµ 22 )).
The PGF of the number of infections of each type emanating from a later-generation infective wearing a mask, given by Γ 1 (s, t), is
where G is the PGF for the excess degree distribution, i.e.,
With the derived PGFs in hand, the probability of extinction starting from a later-generation infective, given by q 1 (resp., q 2 ) if patient zero wears a mask (resp., does not wear a mask), is the smallest non-negative solution of the fixed-point equation (s, t) = (Γ 1 (s, t), Γ 2 (s, t)). Finally, the probability of emergence starting from patient zero, denoted by P 1 (resp., P 2 ) if patient zero wears a mask (resp., does not wear a mask) is given by (P 1 , P 2 ) = (γ 1 (q 1 , q 2 ), γ 2 (q 1 , q 2 )).
Since the probability of emergence is the same in the Mask model and the multi-strain model with mutation, the critical threshold at which an epidemic emerges (also known as the reproductive number R 0 ) is the same in both models as well. To calculate this threshold, we first introduce two matrices:
Then the formula for the critical threshold [8] , [2] is
where ρ(Qµ) denotes the spectral radius of Qµ and k , k 2 are the first and second moments of the degree distribution, respectively. If R 0 < 1 the epidemic dies out in finite time, and if R 0 > 1 an epidemic persists. To write things in terms of the parameters of the Mask model, we can define the matrices
and note that Qµ = Tm. Hence we equivalently have
We follow the method of Eletreby et al [8] . Since the contact network G is drawn from the configuration model with degree distribution {p k } k , it is locally tree-like. We can compute the probability that a given vertex is infected by considering the tree-like neighborhood around it. Mathematically, we can consider an infinite rooted tree where the bottom level is labeled level zero and the top (the root) is labeled level infinity. We let q ,1 (respectively, q ,2 ) be the probability that a mask-wearing (respectively, non-maskwearing) vertex is infected in level .
The pair (q +1,1 , q +1,2 ) can be recursively computed from (q ,1 , q ,2 ) as follows. Consider a vertex in level + 1 that wears a mask. It has degree k with probability kp k k due to properties of the configuration model. Due to the tree structure, k − 1 of these edges are sent to and one is sent to the parent in level + 2. Out of the k − 1 levelneighbors, some number X wear a mask while the rest do not, where X ∼ Binomial(k − 1, m). Out of the X maskwearing neighbors, some number U are infected, where U ∼ Binomial(X, q ,1 ). Similarly, there are V infected non-maskwearing neighbors, where V ∼ Binomial(k − 1 − X, q ,2 ). Finally, if there are U infected mask-wearing neighbors and V infected non-mask-wearing neighbors, the probability that the mask-wearing parent in level + 1 becomes infected is
then we have
If we define f 2 to be the same as (8) except the term
is replaced by
then we also have
Following the analysis in [8] , the sequence {q ,1 , q ,2 } ≥1 converges to a limit (q ∞,1 , q ∞,2 ) which satisfies the fixed point equation
Finally, to compute the probability of infection at the root, we note that the root has k neighbors with probability p k , and all neighbors are in a lower level. Thus, if S 1 (S 2 ) is the probability of infection of a mask-wearing (non-maskwearing) root vertex, then we have
As our analysis was for an arbitrarily chosen root node, S i is the expected fraction of mask-wearing vertices that eventually get infected by the epidemic, conditioned on the epidemic occurring. The total fraction of infections is then given by S = S 1 · m + S 2 · (1 − m).
A. Epidemic as a function of the mean degree We conducted extensive numerical simulations to validate our theoretical analysis. In Figure 1 , we study the probability of emergence. The contact network was generated via the configuration model with Poisson degree distribution and 500,000 vertices. We studied several values for the mean degree ranging between 0 and 10. To generate the simulation plots, we took an average over 20,000 independent trials where, in each trial, a new contact network was generated. The parameters of the mask model were chosen to be m = 0.45, T 11 = 0.126, T 12 = 0.18, T 21 = 0.42, T 22 = 0.6. The choice of m was based on the current fraction of maskwearers in the US [1] . The transmissibility parameters were chosen as a reasonable baseline to illustrate the model and our theoretical results about the model. For larger mean degrees, we see that we have a near-perfect match between the simulations and theoretical predictions. For smaller mean degrees, the match is close, but not perfect, since the emergence event becomes quite rare close to the phase transition point. We expect that if much larger networks are used, the simulations will enjoy better alignment with the theoretical predictions, even close to the phase transition point. In Figure 2 , we study the expected size of the epidemic, conditioned on emergence. In our simulations, we used the same number of nodes and degree distribution, averaged over 10,000 independent trials. The same parameters for the mask model were used as well. We see very good alignment between the simulations and theoretical predictions, confirming the validity of our theoretical results.
In Figure 3 , we illustrate the interesting finding that while the multi-strain model with mutations can be used to compute the probability of emergence in the Mask model, it yields an incorrect prediction of the expected epidemic size. There seems to be a good alignment between the two curves close to the critical threshold, but the two predictions diverge for larger mean degrees. We give a possible reason for this mismatch. In the mask model, there is a single strain in the population and a susceptible vertex is infected as long as as there is a successful infection by at least one neighbor. In the multi-strain model with mutation, if there are multiple successful infections to a susceptible vertex, the resulting transmitted strain depends on the number of successful infections of each type. When the mean degree Fig. 3 . Comparison of the theoretical prediction for the expected epidemic size conditioned on emergence in the mask model (circles) and the multistrain model with mutation (dashed lines). There is good alignment close to the critical threshold, but the predictions diverge for larger mean degrees. We use the same parameters as in Figure 2. is small, it is unlikely that there will be more than one successful infection, as the number of neighbors of a vertex is small. However as the mean degree increases, the difference becomes more pronounced. We plan to further investigate the fundamental differences between the mask model and multi-strain model with mutations in future work. Figure 4 illustrates the effect of the probability of maskwearing, m, on the expected epidemic size. In our simulations, we generated the contact network with 5,000,000 vertices and a Poisson(5) degree distribution. We studied various values of m between 0 and 1. As the fraction of mask-wearing individuals increases, the total number of infections (shown in red) is monotonically decreasing, demonstrating the effectiveness of masks in curbing the spread of COVID-19. Interestingly, we see that the fraction of infected non-mask-wearers (shown in blue) is also monotonically decreasing in m. The intuition for this observation is clear; if many individuals wear a mask, on a high level it reduces the effective transmissibility of the virus, thus reducing the number of infected non-mask-wearers as well. Curiously, the fraction of infected mask-wearers is not monotonically decreasing in m; the infection curve peaks at m = 0.6. We provide a possible explanation. There are two opposing effects which influence the number of infected mask-wearers. As m increases, the total number of infected mask-wearers will naturally increase, since there are more susceptible mask-wearers in the population. On the other hand, increasing m will also decrease the transmissibility of the virus, leading to a lower rate of infection. When m < 0.6, the first effect dominates: the increase in susceptible maskwearers is greater than the decrease in transmissibility. The point m = 0.6 is where the two effect balance each other; for m > 0.6, the decrease in transmissibility dominates the increase in susceptible mask-wearers. Fig. 5 . Plot of the simulated and theoretical critical degree for an epidemic to persist. In our simulations, we generated a contact network with Poisson degree distribution on 5,000 vertices. We averaged over 1,000 independent simulations to generate the simulation data points for the plot.
In Figure 5 , we study how the critical threshold depends on m. As one may expect, as the fraction of mask-wearers increase, a larger mean degree is required for an epidemic to emerge.
In Figure 6 , we conside the effect of the baseline transmissibility (i.e., the transmissibility between two non-maskwearers) on the probability of emergence and expected epidemic size. Instead of setting specific values for the transmissibilities in the mask model, we assume that masks have an inward efficiency of T mask,1 and an outward efficiency of T mask,2 . This implies that the transmission parameters have the form T 11 = T mask,1 T mask,2 T , T 12 = T mask,2 T and T 22 = T . Here, we fix T mask,1 = 0.3 and T mask,2 = 0.7 (this reflects the observation that masks have higher outward Empirical and theoretical plots for the probability of emergence (a) and expected epidemic size (b) as a function of T . The green (blue) curves assume that patient zero is wearing (not wearing) a mask. The red curve assumes that patient zero is randomly selected. efficiency than inward) and we study how the epidemic characteristics change with T . In our simulations, we set m = 0.45 and assumed a Poisson(5) degree distribution and generated networks with 5,000,000 vertices, averaging over 100 independent simulations. In both the probability of emergence and the expected epidemic size, the curves are increasing with T , and an epidemic emerges when T = 0.3. While the simulated probability of emergence deviates from the theoretical curve, we expect to see concentration as we increase the number of experiments.
In this paper, we studied the effects of mask-wearing on viral spread, specifically the probability of emergence and the expected epidemic size conditioned on emergence. We offered two different perspectives on modeling viral spread with masks: through a heterogeneous bond percolation approach on multi-type networks and through an analogy with a multiple-strain model with mutation. Theoretically, we find that while the probability of emergence is the same in both models, the expected epidemic size can be different. We also show that the expected epidemic size is decreasing as a function of the fraction of mask-wearing individuals, confirming that mask-wearing can be an effective strategy in curbing the spread of COVID-19.