# A general epidemic model and its application to mask design considering
different preferences towards masks

Chaoqian Wang and Hamdi Kavak

###### Abstract

While most masks have a limited effect on personal protection, how effective
are they for collective protection? How to enlighten the design of masks from
the perspective of collective dynamics? In this paper, we assume three
preferences in the population: (i) never wearing a mask; (ii) wearing a mask
if and only if infected; (iii) always wearing a mask. We study the epidemic
transmission in an open system within the Susceptible-Infected-Recovered (SIR)
model framework. We use agent-based Monte Carlo simulation and mean-field
differential equations to investigate the model, respectively. Ternary heat
maps show that wearing masks is always beneficial in curbing the spread of the
epidemic. Based on the model, we investigate the potential implications of
different mask designs from the perspective of collective dynamics. The
results show that strengthening the filterability of the mask from the face to
the outside is more effective in most parameter spaces, because it acts on
individuals with both preferences (ii) and (iii). However, when the fraction
of individuals always wearing a mask achieves a critical point, strengthening
the filterability from outside to the face becomes more effective, because of
the emerging hidden reality that the infected individuals become too few to
utilize the filterability from their face to outside fully.

Keywords: Epidemic model; Mask; COVID-19; Verification and validation;
Lyapunov function

##  1 Introduction

As the COVID-19 pandemic is ravaging the world, the protection of masks is a
topic of interest. A news article published in Nature indicates that wearing a
surgical mask leads to an $11\%$ drop in risk, while a $5\%$ drop for cloth
[1]. The protective effect of masks on individuals may seem minimal, but it is
also necessary to focus on the protective effect on collectives.

Since Kermack and McKendrick [2] proposed the Susceptible-Infected-Recovered
(SIR) compartment model, various epidemic models have been developed
considerably. In the classic SIR model, the population is divided into three
compartments: (i) the susceptible ($S$); (ii) the infected ($I$); (ii) the
recovered ($R$). Through human-to-human contact or self-healing, individuals
flow from one compartment to another. A simple modified version is the SEIR
model, which adds an exposed ($E$) compartment to the SIR model. Recently,
Barlow et al. [3, 4] derived the analytical solutions of the SIR [3] and the
SEIR [4] models.

Researchers, in recent years, have explored additional factors and mechanisms
to the classic epidemic models, such as isolation [5] and vaccination [6, 7,
8, 9, 10, 11]. The dynamics of the epidemic transmission can also be applied
to the information spreading, creating rumor spreading models [12, 13] or the
public opinion dynamics model [14, 15]. From the perspective of verification
and validation, the global stability of this class of nonlinear dynamical
systems is widely studied [16, 17, 18, 19, 20, 21]. In particular, researchers
have proved the global stability of endemic equilibria in various epidemic
models in multigroup populations [19, 20, 21], which are general cases of the
model proposed in this work. A common approach to prove global stability is
constructing a Lyapunov function (not limited to epidemiology, but also widely
applied to other complex systems such as evolutionary dynamics [22, 23]),
which measures the system’s “energy.” If the energy continues to decay, then
the system will stabilize at an equilibrium point.

When it comes to the protective effect of masks, several works [24, 25, 26,
27, 28] are noticed to have emerged in the COVID-19 period after 2020. Li et
al. [24] treated whether people wear masks or not as an evolutionary game.
Gondim [25] considered masks in the SEIR model and validated the model by
real-world data. Auger and Moussaoui [26] studied the confinement’s release
threshold, taking the masks into account. Lasisi and Adeyemo [27] modeled the
effect of wearing masks on COVID-19 infection dynamics. Han et al. [28]
investigated the effect of three different preferences on wearing a mask.

Based on the existing literature, we find the previous works on masks have
three shortcomings. First, when classifying the population into three
categories with different preferences on wearing masks according to their
assumptions, there is no work classifying them into three independent
variables. They set only two variables as the fraction of two categories, and
the remaining category’s fraction is one minus these two variables. This leads
to an inability to ensure constant relative proportions of the other two
categories when investigating the effect of the proportion of a certain
category. Second, previous work did not carry out a complete analysis of the
stability of their models. This makes verification and validation challenging.
Third, only focusing on the effect of masks on epidemic spreading, there is no
previous work considering providing applications of the epidemic models to the
design of masks itself.

This paper builds a general epidemic model in an open system considering three
different preferences on wearing masks. We start from a set of agent-based
rules, and use mean-field analysis to verify and validate the model. In
addition to filling in the gaps of previous work by treating three preferences
as independent variables and considering global stability analysis, we explore
[29, 30] the effect of different preferences towards wearing masks on the
epidemic transmission through our model’s eyes. Considering that in the
traditional perception, masks are designed at an individual level, we also try
to reveal the design strategies of the masks by the collective dynamics based
on our model.

##  2 Model

There is an epidemic disease spreading in the system. To prevent this
epidemic, individuals hold different preferences for wearing masks. Concerning
the infection state, we divide the population into: (i) the susceptible ($x$);
(ii) the infected ($y$); (iii) the recovered ($z$). In terms of different
preferences towards wearing masks, we divide the population into: (i) those
who never wear masks (subscript 0); (ii) those who wear masks if and only if
infected (subscript 1); (iii) those who always wear masks (subscript 2).
Therefore, we have up to 9 categories according to different combinations of
the classification of the two dimensions mentioned above.

Before describing evolutionary rules, we list the definition of our
mathematical symbols in Table 1.

Table 1: The definition of mathematical symbols Symbol | Definition
---|---
$x_{0}$ | The number of susceptible individuals never wearing a mask.
$y_{0}$ | The number of infected individuals never wearing a mask.
$z_{0}$ | The number of recovered individuals never wearing a mask.
$x_{1}$ | The number of susceptible individuals wearing a mask if and only if infected.
$y_{1}$ | The number of infected individuals wearing a mask if and only if infected.
$z_{1}$ | The number of recovered individuals wearing a mask if and only if infected.
$x_{2}$ | The number of susceptible individuals always wearing a mask.
$y_{2}$ | The number of infected individuals always wearing a mask.
$z_{2}$ | The number of recovered individuals always wearing a mask.
$n$ | The number of individuals in the system.
$\Lambda$ | The number of new individuals entering the system within unit time.
$\mu$ | The rate of natural death.
$r$ | The rate of recovering.
$\alpha$ | The rate of human-to-human infection.
$\varepsilon_{0}$ | The fraction of new individuals never wearing a mask.
$\varepsilon_{1}$ | The fraction of new individuals wearing a mask if and only if infected.
$\varepsilon_{2}$ | The fraction of new individuals always wearing a mask.
$p_{I}$ | The protective effect produced when an infected individual wears a mask.
$p_{S}$ | The protective effect produced when a susceptible individual wears a mask.

###  2.1 The agent-based rules

Consider an open system containing initially $n\big{|}_{t=0}$ agents (i.e.,
individuals). Within a Monte Carlo step, an agent $i$ is randomly selected,
and the following parallel events occur.

(1) If agent $i$ is susceptible, we again select an agent $j$ randomly. If
agent $j$ is infected, then agent $i$ is infected with a probability $\alpha$
($\alpha>0$). If agent $j$ wears a mask, then agent $i$ spares from infection
with a probability $p_{I}$ ($0<p_{I}<1$). If agent $i$ wears a mask, then
agent $i$ spares from infection with a probability $p_{S}$ ($0<p_{S}<1$).

(2) If agent $i$ is infected, then it recovers with a probability $r$ (the
average infection cycle is $1/r$). This does not happen at the same Monte
Carlo step as the event (1).

(3) Agent $i$ naturally dies with a probability $\mu$ (the average lifespan is
$1/\mu$). We do not consider deaths due to the epidemic.

To ensure the population remains almost unchanged, we must let new agents
enter the system. We set the following very first event in a Monte Carlo step,
where the number (0) means it happens before the event (1).

(0) A new agent enters the system with a probability $p_{e}$. The agent’s
personal preference determines it never wears a mask with a probability
$\varepsilon_{0}$ ($0<\varepsilon_{0}<1$), wears a mask if and only if
infected with a probability $\varepsilon_{1}$ ($0<\varepsilon_{1}<1$), or
always wears a mask with a probability $\varepsilon_{2}$
($0<\varepsilon_{2}<1$), yielding
$\varepsilon_{0}+\varepsilon_{1}+\varepsilon_{2}=1$. The preference of an
agent on masks does not change over time.

For the population $n$ remaining almost unchanged with time, we let one time
step $t$ contains $n\big{|}_{t=0}$ Monte Carlo steps, such that each agent can
be selected once on average. Therefore, the expected number of new agents
entering the system within one time step $t$ is $\Lambda=n^{*}p_{e}$. The
solution is $p_{e}=\mu$ (see Theorem 1).

###  2.2 The mean-field equations

Performing mean-field analysis, we can approximate the agent-based dynamics
into a set of differential equations. We do not dwell on the principles of
mean-field analysis, but only explain some important points. (i) Within unit
time, each agent is selected once on average, such that the number of events
descending in a category is the population in the category. (ii) The
probability of selecting an infected agent never wearing a mask from the
population is $y_{0}/n$, and $y_{1}/n$, $y_{2}/n$ for selecting an infected
agent holding the other two preferences, respectively. (iii) Thanks to the
mask, the probability of sparing from infection is $p_{I}$ or $p_{S}$, which
means the probability of infection is $(1-p_{I})$ or $(1-p_{S})$. (iv) If
there are two layers of protection, they must be both breached for the
infection to succeed.

We denote the state of the system by a vector $\mathbf{\Psi}$, containing the
population in nine categories. The mean-field differential equations depicting
the agent-based dynamics is

$\dot{\mathbf{\Psi}}=\begin{pmatrix}\dot{x}_{0}\\\ \dot{y}_{0}\\\
\dot{z}_{0}\\\ \dot{x}_{1}\\\ \dot{y}_{1}\\\ \dot{z}_{1}\\\ \dot{x}_{2}\\\
\dot{y}_{2}\\\ \dot{z}_{2}\end{pmatrix},$ (1)

where

$\left\\{\begin{aligned} \dot{x}_{0}=&~{}\varepsilon_{0}\Lambda-\alpha
x_{0}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n-\mu x_{0},\\\ \dot{y}_{0}=&~{}\alpha
x_{0}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n-ry_{0}-\mu y_{0},\\\
\dot{z}_{0}=&~{}ry_{0}-\mu z_{0},\\\
\dot{x}_{1}=&~{}\varepsilon_{1}\Lambda-\alpha
x_{1}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n-\mu x_{1},\\\ \dot{y}_{1}=&~{}\alpha
x_{1}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n-ry_{1}-\mu y_{1},\\\
\dot{z}_{1}=&~{}ry_{1}-\mu z_{1},\\\
\dot{x}_{2}=&~{}\varepsilon_{2}\Lambda-\alpha(1-p_{S})x_{2}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n\\\
&-\mu x_{2},\\\
\dot{y}_{2}=&~{}\alpha(1-p_{S})x_{2}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n-ry_{2}\\\
&-\mu y_{2},\\\ \dot{z}_{2}=&~{}ry_{2}-\mu z_{2}.\end{aligned}\right.$

The system depicted by Eq. (1) is an extended version of the SIR model with a
standard incidence rate.

##  3 Results and discussion

In this section, we demonstrate the numerical results of the model from both
Monte Carlo simulation and mean-field equations. The algorithm of the Monte
Carlo simulation was described in Section 2.1. In the numerical simulation of
mean-field equations, we use iteration $\mathbf{\Psi}(t+\Delta
t)=\mathbf{\Psi}(t)+\dot{\mathbf{\Psi}}(t)\Delta t$, where $\Delta t=0.01$.

We set three statistical measures: (i) the proportion of susceptible,
$p_{x}=(x_{0}+x_{1}+x_{2})/n$; (ii) the proportion of infected,
$p_{y}=(y_{0}+y_{1}+y_{2})/n$; (iii) the proportion of recovered,
$p_{z}=(z_{0}+z_{1}+z_{2})/n$.

###  3.1 Impact of masks on epidemic spreading

Figures 1 and 2 show the time evolution of $p_{x}$, $p_{y}$, and $p_{z}$ with
different parameters and initial conditions.

Figure 1:  Time evolution of $p_{x}$, $p_{y}$, and $p_{z}$. The results of
Monte Carlo simulation and mean-field equations are presented together.
$\alpha=0.2$, $\mu=0.01$, $\Lambda=10$, $r=0.05$, $\varepsilon_{0}=0.1$,
$\varepsilon_{1}=0.1$, $\varepsilon_{2}=0.8$, $p_{I}=0.75$, $p_{S}=0.25$. (a)
$p_{x}\big{|}_{t=0}=0.9$, $p_{y}\big{|}_{t=0}=0.1$, $p_{z}\big{|}_{t=0}=0$.
(b) $p_{x}\big{|}_{t=0}=0.1$, $p_{y}\big{|}_{t=0}=0.9$,
$p_{z}\big{|}_{t=0}=0$. (a)(b)
$x_{0}\big{|}_{t=0}=x_{1}\big{|}_{t=0}=x_{2}\big{|}_{t=0}$ in
$x\big{|}_{t=0}$, as well as $y\big{|}_{t=0}$ and $z\big{|}_{t=0}$. Figure 2:
Time evolution of $p_{x}$, $p_{y}$, and $p_{z}$. The results of Monte Carlo
simulation and mean-field equations are presented together. $\alpha=0.2$,
$\mu=0.01$, $\Lambda=10$, $r=0.05$, $\varepsilon_{0}=0.3$,
$\varepsilon_{1}=0.1$, $\varepsilon_{2}=0.6$, $p_{I}=0.5$, $p_{S}=0.05$. (a)
$p_{x}\big{|}_{t=0}=0.9$, $p_{y}\big{|}_{t=0}=0.1$, $p_{z}\big{|}_{t=0}=0$.
(b) $p_{x}\big{|}_{t=0}=0.1$, $p_{y}\big{|}_{t=0}=0.9$,
$p_{z}\big{|}_{t=0}=0$. (a)(b)
$x_{0}\big{|}_{t=0}=x_{1}\big{|}_{t=0}=x_{2}\big{|}_{t=0}$ in
$x\big{|}_{t=0}$, as well as $y\big{|}_{t=0}$ and $z\big{|}_{t=0}$.

From Figures 1 and 2, we find the proportions of different individuals always
achieve stability after time evolution. The results of the Monte Carlo
simulation fluctuate, while the results of mean-field equations are stable.
They corroborate each other. In addition, we obverse two phenomena. First, the
epidemic may either die out or exist at the end, dependent on different
parameters. Second, with the same parameters and different initial conditions,
the steady-states are the same.

Next, in the heat maps Figure 3 and 4, we present the steady-states of $p_{y}$
and $p_{x}$ as a ternary function of $\varepsilon_{0}$, $\varepsilon_{1}$,
$\varepsilon_{2}$, respectively. The results of the Monte Carlo simulation are
the average of the last 200 time steps ($t$). The results of mean-field
equations are retrieved from the state $\mathbf{\Psi}(t+\Delta t)$ when
$\max\\{|\mathbf{\Psi}(t)|\Delta t\\}<0.0001$.

Figure 3:  The steady-state of $p_{y}$ as a ternary function of
$\varepsilon_{0}$, $\varepsilon_{1}$, $\varepsilon_{2}$. (a) Monte Carlo
simulation. (b) Mean-field equations. $\alpha=0.2$, $\mu=0.01$, $\Lambda=10$,
$r=0.05$, $p_{I}=0.5$, $p_{S}=0.05$. Figure 4:  The steady-state of $p_{x}$ as
a ternary function of $\varepsilon_{0}$, $\varepsilon_{1}$, $\varepsilon_{2}$.
(a) Monte Carlo simulation. (b) Mean-field equations. $\alpha=0.2$,
$\mu=0.01$, $\Lambda=10$, $r=0.05$, $p_{I}=0.5$, $p_{S}=0.05$.

From Figures 3 and 4, we find that the results of Monte Carlo simulation and
mean-field equations corroborate each other. In Figure 3, we observe that more
individuals wearing masks reduce the proportion of infected individuals in the
population. In particular, always wearing a mask has a better effect on
reducing infected individuals. In Figure 4, we observe that more individuals
wearing masks increases the proportion of susceptible individuals in the
population. Different from increasing recovered individuals, it means that
more people are spared from getting infected once. Also, always wearing a mask
has a better effect on increasing susceptible individuals (Figure 4), making
more individuals spared from being infected.

###  3.2 Potential implication of different mask design

Based on our model, we can reveal the potential implication of different mask
designs. When an infected individual wears a mask, it is the filterability
from the face to the outside that provides protection (to the population), and
when a susceptible individual wears a mask, the filterability from the outside
to the face matters.

We show in Figure 5 (a) and (b) the steady infected fraction $p_{y}$ as a
binary function of the protective effect produced when an infected individual
wears a mask ($p_{I}$) and when a susceptible individual wears a mask
($p_{S}$).

Figure 5:  (a) Monte Carlo simulation. The steady-state of $p_{y}$ as a binary
function of $p_{I}$, $p_{S}$. (b) Mean-field equations. The steady-state of
$p_{y}$ as a binary function of $p_{I}$, $p_{S}$. (c) Mean-field equations.
The steady-state of $\partial p_{y}/\partial p_{I}$ and $\partial
p_{y}/\partial p_{S}$ as a binary function of $p_{I}$, $p_{S}$. $\alpha=0.2$,
$\mu=0.01$, $\Lambda=10$, $r=0.05$, $\varepsilon_{0}=0.3$,
$\varepsilon_{1}=0.1$, $\varepsilon_{2}=0.6$, $p_{I}=0.5$, $p_{S}=0.05$.

We can observe that in both the Monte Carlo simulation [Figure 5(a)] and mean-
field equations [Figure 5(b)], an increase in protective effect $p_{I}$ or
$p_{S}$ leads to a decrease in the steady infected fraction $p_{y}$. On this
basis, we further ask which one in increasing $p_{I}$ or $p_{S}$ is more
effective in reducing the infected fraction?

We show in Figure 5(c) the steady-state of $\partial p_{y}/\partial p_{I}$ and
$\partial p_{y}/\partial p_{S}$ as a binary function of $p_{I}$ and $p_{S}$.
If $\partial p_{y}/\partial p_{I}<\partial p_{y}/\partial p_{S}$, increasing
the unit protective effect from the infected side is more conducive to
reducing the infected fraction, and vice versa. Intuitively, increasing
$p_{I}$ should have always been more conducive than increasing $p_{S}$,
because the former acts on individuals with two preferences, (ii) those who
wear masks if and only if infected and (iii) those who always wear masks. In
contrast, the latter only acts on individuals with one preference, (iii) those
who always wear masks. Increasing $p_{I}$ obviously has a broader scope of
action than the increasing $p_{S}$ and covers the latter’s population.
However, Figure 5(c) presents a different phenomenon. This indicates that we
can provide the designs of the masks with different insights from the group
dynamics. We will give this further analysis in Section 4.5.

##  4 Verification and validation

This section verifies and validates the properties that we concluded in
numerical results by analyzing them at a mathematical level.

###  4.1 The total population dynamics

The total population dynamics follows Theorem 1.

######  Theorem 1.

For $t\to\infty$, we have $n\to\Lambda/\mu$.

######  Proof.

According to Eq. (1),

$\displaystyle\dot{n}$
$\displaystyle=\dot{x_{0}}+\dot{y_{0}}+\dot{z_{0}}+\dot{x_{1}}+\dot{y_{1}}+\dot{z_{1}}+\dot{x_{2}}+\dot{y_{2}}+\dot{z_{2}}$
$\displaystyle=\Lambda-\mu n.$ (2)

Solving Eq. (Proof), we get

$n=\left(n\big{|}_{t=0}-\frac{\Lambda}{\mu}\right)\mathrm{e}^{-\mu
t}+\frac{\Lambda}{\mu}.$ (3)

From Eq. (3), we complete the proof that $n\to\Lambda/\mu$ for $t\to\infty$.
We give this significant value a symbol $n^{*}$,

$n^{*}=\lim_{t\to\infty}n=\frac{\Lambda}{\mu}.$ (4)

In the same way, we can also prove that
$x_{0}+y_{0}+z_{0}\to\varepsilon_{0}\Lambda/\mu$,
$x_{1}+y_{1}+z_{1}\to\varepsilon_{1}\Lambda/\mu$,
$x_{2}+y_{2}+z_{2}\to\varepsilon_{2}\Lambda/\mu$ for $t\to\infty$.

Theorem 1 gives us another important insight: when discussing the steady-
state, we can substitute $n$ for $n^{*}=\Lambda/\mu$ in the system of Eq. (1).

###  4.2 The basic reproduction number

The basic reproduction number $\mathcal{R}_{0}$ is one of the most important
measures in epidemiology. It can assist in analyzing both the stability of the
system and the effect of parameters on the epidemic transmission.

We let $\dot{\mathbf{\Psi}}=\mathbf{0}$ and solve for the epidemic-free
equilibrium, denoted by $\mathbf{\Psi}^{*}$,

$\mathbf{\Psi}^{*}=\frac{\Lambda}{\mu}\left(\varepsilon_{0},0,0,\varepsilon_{1},0,0,\varepsilon_{2},0,0\right)^{\mathrm{T}}.$
(5)

Then, we can follow the method in Ref. [31] to find the basic reproduction
number (see Appendix A):

$\mathcal{R}_{0}=\frac{\alpha}{r+\mu}[\varepsilon_{0}+(1-p_{I})\varepsilon_{1}+(1-p_{S})(1-p_{I})\varepsilon_{2}].$
(6)

The basic reproduction number reveals the following theorem.

######  Theorem 2.

The epidemic-free equilibrium $\mathbf{\Psi}^{*}$ is locally asymptotically
stable if $\mathcal{R}_{0}<1$, and the epidemic-free equilibrium
$\mathbf{\Psi}^{*}$ is not stable if $\mathcal{R}_{0}>1$.

(See Ref. [31] for proof)

Substituting the parameters in Figure 1 into Eq. (6), we can calculate
$\mathcal{R}_{0}=0.9167<1$, which means the epidemic-free equilibrium is
stable, consistent with that shown in Figure 1. Similarly, substituting
parameters in Figure 2 produces $\mathcal{R}_{0}=2.1167>1$, such that the
epidemic-free equilibrium is not stable, which is also consistent with that
shown in Figure 2.

###  4.3 Epidemic-free and endemic equilibria

We separate $\mathbf{\Psi}$ into

$\left\\{\begin{aligned}
\mathbf{\Phi}_{1}&=\left(x_{0},y_{0},x_{1},y_{1},x_{2},y_{2}\right)^{\mathrm{T}},\\\
\mathbf{\Phi}_{2}&=\left(z_{0},z_{1},z_{2}\right)^{\mathrm{T}}.\end{aligned}\right.$
(7)

From Eq. (1), we can assert that the steady-state of $\mathbf{\Phi}_{1}$ can
determine the steady-state of $\mathbf{\Phi}_{2}$, and $\mathbf{\Phi}_{2}$
does not affect the evolution of $\mathbf{\Phi}_{1}$. Therefore, the stability
of $\mathbf{\Psi}$ is equivalent to (i) the stability of $\mathbf{\Phi}_{1}$
and (ii) the stability of $\mathbf{\Phi}_{2}$ when $\mathbf{\Phi}_{1}$
achieves stability.

We first prove the global stability of the epidemic-free equilibrium by
constructing the Lyapunov function.

######  Theorem 3.

The epidemic-free equilibrium $\mathbf{\Psi}^{*}$ is global asymptotically
stable in $\mathbb{R}_{\geq 0}^{9}$ if $\mathcal{R}_{0}<1$.

######  Proof.

Consider the Lyapunov function $\mathcal{L}(\mathbf{\Phi}_{1})$ in
$\mathbb{R}_{\geq 0}^{6}$,

$\displaystyle\mathcal{L}(\mathbf{\Phi}_{1})=$
$\displaystyle~{}\frac{(x_{0}-x_{0}^{*})^{2}}{2x_{0}^{*}}+y_{0}$
$\displaystyle+(1-p_{I})\left[\frac{(x_{1}-x_{1}^{*})^{2}}{2x_{1}^{*}}+y_{1}\right]$
$\displaystyle+(1-p_{I})\left[\frac{(x_{2}-x_{2}^{*})^{2}}{2x_{2}^{*}}+y_{2}\right].$
(8)

We can conclude that, (i) $\mathcal{L}(\mathbf{\Phi}_{1})=0$ when
$\mathbf{\Phi}_{1}=\mathbf{\Phi}_{1}^{*}$, (ii)
$\mathcal{L}(\mathbf{\Phi}_{1})>0$ when
$\mathbf{\Phi}_{1}\neq\mathbf{\Phi}_{1}^{*}$. Therefore, $\mathbf{\Phi}_{1}$
is positive definite in the neighborhood of $\mathbf{\Phi}_{1}^{*}$. Secondly,
we have (see Appendix B),

$\dot{\mathcal{L}}(\mathbf{\Phi}_{1})\leq 0$ (9)

when $\mathbf{\Phi}_{1}\neq\mathbf{\Phi}_{1}^{*}$. Note that
$\dot{\mathcal{L}}(\mathbf{\Phi}_{1})=0$ can be confirmed by Eq. (B) when
$\mathbf{\Phi}_{1}=\mathbf{\Phi}_{1}^{*}$. Therefore, $\mathbf{\Phi}_{1}$ is
negative semi-definite in the neighborhood of $\mathbf{\Phi}_{1}^{*}$. Hence,
according to Lasalle’s Invariance Principle [32], $\mathbf{\Phi}_{1}^{*}$ is
globally asymptotically stable in $\mathbb{R}_{\geq 0}^{6}$. Given
$\mathbf{\Phi}_{1}^{*}$ stable, it is easy to prove the global asymptotic
stability of $\mathbf{\Phi}_{2}^{*}$ in $\mathbb{R}_{\geq 0}^{3}$ by
constructing Lyapunov function
$\mathcal{L}_{z}(\mathbf{\Phi}_{2})=z_{0}+z_{1}+z_{2}$. Therefore,
$\mathbf{\Psi}^{*}$ is global asymptotically stable in $\mathbb{R}_{\geq
0}^{9}$.

The equation $\dot{\mathbf{\Psi}}=\mathbf{0}$ has two solutions. We denote the
second solution by $\mathbf{\Psi}^{**}$. The infected population is non-zero;
thus, we call it the endemic equilibrium. This equilibrium
$\mathbf{\Psi}^{**}$ corresponds to the results shown in Figure 2. It is not
easy to express it analytically, but we can still show its existence
condition.

######  Theorem 4.

The endemic equilibrium $\mathbf{\Psi}^{**}$ exists and is unique in
$\mathbb{R}_{\geq 0}^{9}$ if $\mathcal{R}_{0}>1$.

######  Proof.

First, we show the relationship between the existence and uniqueness of
positive $y_{i}^{**}$, $i=0,1,2$ ($y_{0}^{**}>0$, $y_{1}^{**}>0$,
$y_{2}^{**}>0$) and $\mathcal{R}_{0}>1$ (see Appendix C). Then, the existence
and uniqueness of $x_{i}^{**}$, $z_{i}^{**}$, $i=0,1,2$ can be naturally
confirmed, hence the existence and uniqueness of $\mathbf{\Psi}^{**}$.

Theorem 3 validates that in Figure 1, the steady-state with the same
parameters is independent of the initial conditions. Theorem 4 is consistent
with Figure 2.

###  4.4 Robustness analysis for the effect of wearing masks

The basic reproduction number measures the average number of individuals that
an infected individual can transmit the epidemic. The higher the basic
reproduction number, the more severe the epidemic.

Analyzing Eq. (6), we can see that the coefficient before $\varepsilon_{0}$ is
1, and $(1-p_{I})$ for $\varepsilon_{1}$, and $(1-p_{S})(1-p_{I})$ for
$\varepsilon_{2}$. Since $p_{S}>0$, $p_{I}>0$, we have
$(1-p_{S})(1-p_{I})<1-p_{I}<1$. Considering the constraints:
$0\leq\varepsilon_{0}\leq 1$, $0\leq\varepsilon_{1}\leq 1$,
$0\leq\varepsilon_{2}\leq 1$,
$\varepsilon_{0}+\varepsilon_{1}+\varepsilon_{2}=1$, we know the following
facts. (i) $\mathcal{R}_{0}$ takes the minimum when $\varepsilon_{0}=0$,
$\varepsilon_{1}=0$, $\varepsilon_{2}=1$. (ii) $\mathcal{R}_{0}$ takes the
maximum when $\varepsilon_{0}=1$, $\varepsilon_{1}=0$, $\varepsilon_{2}=0$.
Therefore, everyone always wearing a mask minimizes the epidemic severity,
while no one wearing masks maximizes the epidemic severity.

###  4.5 Application to mask design

The basic reproduction number $\mathcal{R}_{0}$ can play the same role as the
infected fraction $p_{y}$ in measuring the outbreak severity. From Eq. (6), we
see that the coefficient $(1-p_{I})$ acts on both $\varepsilon_{1}$ and
$\varepsilon_{2}$ while $(1-p_{S})$ only acts on $\varepsilon_{2}$, which
means the mask protection from the infected side acts on two population
categories and that from the susceptible side acts on only one. This brings us
a misleading intuition that increasing the protective effect from the infected
side is always more conducive. However, we can reveal a hidden different
reality by group dynamics.

We write the partial derivatives of $\mathcal{R}_{0}$ with respect to $p_{I}$
and $p_{S}$ in Eq. (10) and (11).

$\frac{\partial\mathcal{R}_{0}}{\partial
p_{I}}=-\frac{\alpha}{r+\mu}[\varepsilon_{1}+(1-p_{S})\varepsilon_{2}],$ (10)
$\frac{\partial\mathcal{R}_{0}}{\partial
p_{S}}=-\frac{\alpha}{r+\mu}(1-p_{I})\varepsilon_{2}.$ (11)

Increasing the protective effect from the infected side is better than that of
the susceptible one means

$\frac{\partial\mathcal{R}_{0}}{\partial
p_{I}}<\frac{\partial\mathcal{R}_{0}}{\partial p_{S}}$ (12)

or

$\frac{\varepsilon_{1}}{\varepsilon_{2}}>p_{S}-p_{I},$ (13)

and, increasing the protective effect from the susceptible side is better than
that of the infected one means

$\frac{\partial\mathcal{R}_{0}}{\partial
p_{S}}<\frac{\partial\mathcal{R}_{0}}{\partial p_{I}}$ (14)

or

$\frac{\varepsilon_{1}}{\varepsilon_{2}}<p_{S}-p_{I},$ (15)

We can discuss the parameter space in two cases. First, if $p_{S}<p_{I}$, then
Eq. (12) always holds. Second, if $p_{S}>p_{I}$, then Eq. (14) does not always
hold. This suggests that the “intuitive" phenomena (increasing $p_{I}$ is more
effective) occupy more parameter space, which can be verified by Figure 5(c).

For the latter case, $p_{S}>p_{I}$, we can transform Eq. (15) into
$\varepsilon_{2}>\varepsilon_{1}/(p_{S}-p_{I})$. This indicates that if the
fraction of individuals always wearing a mask ($\varepsilon_{2}$) exceeds a
critical point, $\varepsilon_{1}/(p_{S}-p_{I})$, then, increasing $p_{S}$ is
more effective than increasing $p_{I}$, even if $p_{I}$ can act on both
category $\varepsilon_{1}$ and $\varepsilon_{2}$. We can interpret this result
in daily language by taking into account the fraction of existing infected
individuals. As we showed in Section 4.4, $\mathcal{R}_{0}$ decreases (i.e.,
infected individuals increasing) with an increase in $\varepsilon_{2}$. In
this way, the critical point of $\varepsilon_{2}$ is rational to exist, over
which the infected individuals are too few to exert the protective effect that
the mask produces on their side.

##  5 Conclusion

Although most masks have little to no effect on personal protection [1], we
are still interested in the protective effects of masks on a population. We
proposed a general epidemic model in the classic SIR framework considering
three different preferences towards wearing masks. Some individuals never wear
masks; others wear masks if and only if infected, and some always wear masks.
We started from agent-based rules and used a set of mean-field differential
equations to approximate the model. The results of the two corroborate each
other. In this work, the three preferences are independent of each other.

The first aspect is the effect of masks on epidemics. The ternary heat maps
revealed that wearing masks can reduce the number of infected individuals and
increase the number of susceptible individuals. We provided the global
stability analysis of the results and showed the robustness of the
effectiveness of masks by analyzing the basic reproduction number of the
epidemic. We concluded that wearing masks are beneficial to the control of
epidemics.

The second aspect is the application of the epidemic model to mask design. The
protective effect from the infected side ($p_{I}$) can be understood as the
filterability of the mask from the face to the outside against viruses, while
the protective effect to the susceptible side ($p_{S}$) can be interpreted as
the filterability from the outside to the face. This can be influenced by the
material and design of the mask [33], and we analyzed which side strengthening
would provide better results. We showed that strengthening the infected side
is more effective in most parameter spaces. This is intuitive since
strengthening the infected side acts on two categories of individuals (those
wearing masks only if infected and those always wearing masks), while
strengthening the susceptible side acts on only one category (those always
wearing masks). However, there is a hidden reality from the perspective of
group dynamics. We found that once the fraction of individuals always wearing
masks exceeds a critical point,
$\varepsilon_{2}>\varepsilon_{1}/(p_{S}-p_{I})$, then, strengthening the
susceptible side becomes more effective. This is because the preference of
always wearing masks reduces the infected fraction in the population, so that
the infected individuals are too few to exert the protective effect of masks
produced on their side. In the daily language, both the cases above seem to
make sense. However, noticing the latter case from the group perspective and
further giving the mask design strategies according to parameter spaces are
not straightforward without the help of system dynamics.

Real-world situations may have more complexity and different insights. For
instance, the underlying assumptions—people’s preferences do not change with
time, ignores human subjectivity, which has the potential to reveal more
insights. In fact, people can change their preference on whether to wear masks
by either estimating the epidemic severity (evolutionary games) or being
affected by the propaganda of the effectiveness of masks (opinion dynamics).
In this way, future work may consider time-dependent preferences, and apply
any modified model to mask design.

## Acknowledgement

Publication of this article was funded in part by the George Mason University
Libraries Open Access Publishing Fund.

## Data availability

Data sharing not applicable to this article as no datasets were generated or
analyzed during the current study.

## Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

##  Appendix A Finding the basic reproduction number

Let us decompose the infected compartments in Eq. (1) as
$\left(\dot{y}_{0},\dot{y}_{1},\dot{y}_{2}\right)^{\mathrm{T}}=\mathcal{F}-\mathcal{V}$,
where

$\mathcal{F}=\begin{pmatrix}\mathcal{F}_{y_{0}}\\\ \mathcal{F}_{y_{1}}\\\
\mathcal{F}_{y_{0}}\end{pmatrix}=\begin{pmatrix}\alpha
x_{0}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n^{*}\\\ \alpha
x_{1}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n^{*}\\\
\alpha(1-p_{S})x_{2}[y_{0}+(1-p_{I})(y_{1}+y_{2})]/n^{*}\end{pmatrix},$ (A.1)
$\mathcal{V}=\begin{pmatrix}\mathcal{V}_{y_{0}}\\\ \mathcal{V}_{y_{1}}\\\
\mathcal{V}_{y_{0}}\end{pmatrix}=\begin{pmatrix}ry_{0}+\mu y_{0}\\\ ry_{1}+\mu
y_{1}\\\ ry_{2}+\mu y_{2}\end{pmatrix}.$ (A.2)

Solve for the Jacobian matrix of $\mathcal{F}$ and $\mathcal{V}$ at
$\mathbf{\Psi}^{*}$, denoted by $\mathbf{F}$ and $\mathbf{V}$,

$\displaystyle\mathbf{F}$
$\displaystyle=\begin{pmatrix}\displaystyle\frac{\partial\mathcal{F}_{y_{0}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{F}_{y_{0}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{F}_{y_{0}}}{\partial
y_{2}}\\\\[8.0pt] \displaystyle\frac{\partial\mathcal{F}_{y_{1}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{F}_{y_{1}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{F}_{y_{1}}}{\partial
y_{2}}\\\\[8.0pt] \displaystyle\frac{\partial\mathcal{F}_{y_{2}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{F}_{y_{2}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{F}_{y_{2}}}{\partial
y_{2}}\end{pmatrix}(\mathbf{\Psi}^{*})=\frac{1}{n^{*}}\begin{pmatrix}\displaystyle\alpha
x_{0}&\displaystyle\alpha(1-p_{I})x_{0}&\displaystyle\alpha(1-p_{I})x_{0}\\\
\displaystyle\alpha
x_{1}&\displaystyle\alpha(1-p_{I})x_{1}&\displaystyle\alpha(1-p_{I})x_{1}\\\
\displaystyle\alpha(1-p_{S})x_{2}&\displaystyle\alpha(1-p_{S})(1-p_{I})x_{2}&\displaystyle\alpha(1-p_{S})(1-p_{I})x_{2}\end{pmatrix}$
$\displaystyle=\alpha\begin{pmatrix}\displaystyle
x_{0}&\displaystyle(1-p_{I})x_{0}&\displaystyle(1-p_{I})x_{0}\\\ \displaystyle
x_{1}&\displaystyle(1-p_{I})x_{1}&\displaystyle(1-p_{I})x_{1}\\\
\displaystyle(1-p_{S})x_{2}&\displaystyle(1-p_{S})(1-p_{I})x_{2}&\displaystyle(1-p_{S})(1-p_{I})x_{2}\end{pmatrix},$
(A.3)
$\mathbf{V}=\begin{pmatrix}\displaystyle\frac{\partial\mathcal{V}_{y_{0}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{V}_{y_{0}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{V}_{y_{0}}}{\partial
y_{2}}\\\\[8.0pt] \displaystyle\frac{\partial\mathcal{V}_{y_{1}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{V}_{y_{1}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{V}_{y_{1}}}{\partial
y_{2}}\\\\[8.0pt] \displaystyle\frac{\partial\mathcal{V}_{y_{2}}}{\partial
y_{0}}&\displaystyle\frac{\partial\mathcal{V}_{y_{2}}}{\partial
y_{1}}&\displaystyle\frac{\partial\mathcal{V}_{y_{2}}}{\partial
y_{2}}\end{pmatrix}(\mathbf{\Psi}^{*})=(r+\mu)\begin{pmatrix}\displaystyle
1&\displaystyle 0&\displaystyle 0\\\ \displaystyle 0&\displaystyle
1&\displaystyle 0\\\ \displaystyle 0&\displaystyle 0&\displaystyle
1\end{pmatrix}.$ (A.4)

Then, the spectral radius (i.e., maximum eigenvalue) of
$\mathbf{F}\cdot\mathbf{V}^{-1}$ is the basic reproduction number
$\mathcal{R}_{0}$,

$\mathcal{R}_{0}=\frac{\alpha}{r+\mu}[\varepsilon_{0}+(1-p_{I})\varepsilon_{1}+(1-p_{S})(1-p_{I})\varepsilon_{2}].$
(A.5)

Please see Ref. [31] for more information on how to find the basic
reproduction number.

##  Appendix B Proof of $\dot{\mathcal{L}}(\mathbf{\Phi}_{1})\leq 0$ when
$\mathbf{\Phi}_{1}\neq\mathbf{\Phi}_{1}^{*}$

$\displaystyle\dot{\mathcal{L}}(\mathbf{\Phi}_{1})=$
$\displaystyle\left(\frac{x_{0}}{x_{0}^{*}}-1\right)\dot{x}_{0}+\dot{y}_{0}+(1-p_{I})\left[\left(\frac{x_{1}}{x_{1}^{*}}-1\right)\dot{x}_{1}+\dot{y}_{1}\right]+(1-p_{I})\left[\left(\frac{x_{2}}{x_{2}^{*}}-1\right)\dot{x}_{2}+\dot{y}_{2}\right]$
$\displaystyle=$
$\displaystyle\left(\frac{x_{0}}{x_{0}^{*}}-1\right)\left(\varepsilon_{0}\Lambda-\frac{\alpha
x_{0}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}-\mu x_{0}\right)+\frac{\alpha
x_{0}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}$ $\displaystyle-ry_{0}-\mu
y_{0}+(1-p_{I})\left(\frac{x_{1}}{x_{1}^{*}}-1\right)\left(\varepsilon_{1}\Lambda-\frac{\alpha
x_{1}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}-\mu x_{1}\right)$
$\displaystyle+(1-p_{I})\left(\frac{\alpha
x_{1}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}-ry_{1}-\mu y_{1}\right)$
$\displaystyle+(1-p_{I})\left(\frac{x_{2}}{x_{2}^{*}}-1\right)\left(\varepsilon_{2}\Lambda-\frac{\alpha(1-p_{S})x_{2}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}-\mu
x_{2}\right)$
$\displaystyle+(1-p_{I})\left(\frac{\alpha(1-p_{S})x_{2}[y_{0}+(1-p_{I})(y_{1}+y_{2})]}{n^{*}}-ry_{2}-\mu
y_{2}\right)$ $\displaystyle=$
$\displaystyle-\frac{\mu}{x_{0}^{*}}(x_{0}-x_{0}^{*})^{2}-\frac{\alpha}{x_{0}^{*}n^{*}}[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{0}-x_{0}^{*})^{2}$
$\displaystyle+(r+\mu)\left(\frac{\alpha
x_{0}^{*}}{r+\mu}\times\frac{y_{0}+(1-p_{I})(y_{1}+y_{2})}{n^{*}}-y_{0}\right)$
$\displaystyle-\frac{\mu}{x_{1}^{*}}(1-p_{I})(x_{1}-x_{1}^{*})^{2}-\frac{\alpha}{x_{1}^{*}n^{*}}(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{1}-x_{1}^{*})^{2}$
$\displaystyle+(r+\mu)(1-p_{I})\left(\frac{\alpha
x_{1}^{*}}{r+\mu}\times\frac{y_{0}+(1-p_{I})(y_{1}+y_{2})}{n^{*}}-y_{1}\right)$
$\displaystyle-\frac{\mu}{x_{2}^{*}}(1-p_{I})(x_{2}-x_{2}^{*})^{2}-\frac{\alpha}{x_{2}^{*}n^{*}}(1-p_{S})(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{2}-x_{2}^{*})^{2}$
$\displaystyle+(r+\mu)(1-p_{I})\left(\frac{\alpha(1-p_{S})x_{2}^{*}}{r+\mu}\times\frac{y_{0}+(1-p_{I})(y_{1}+y_{2})}{n^{*}}-y_{2}\right).$
(B.1)

In Eq. (B), we used $x_{0}^{*}=\varepsilon_{0}\Lambda/\mu$,
$x_{1}^{*}=\varepsilon_{1}\Lambda/\mu$,
$x_{2}^{*}=\varepsilon_{2}\Lambda/\mu$.

We can further deflate Eq. (B),

$\displaystyle\dot{\mathcal{L}}(\mathbf{\Phi}_{1})=$
$\displaystyle-\frac{\mu}{x_{0}^{*}}(x_{0}-x_{0}^{*})^{2}-\frac{\alpha}{x_{0}^{*}n^{*}}[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{0}-x_{0}^{*})^{2}$
$\displaystyle-\frac{\mu}{x_{1}^{*}}(1-p_{I})(x_{1}-x_{1}^{*})^{2}-\frac{\alpha}{x_{1}^{*}n^{*}}(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{1}-x_{1}^{*})^{2}$
$\displaystyle-\frac{\mu}{x_{2}^{*}}(1-p_{I})(x_{2}-x_{2}^{*})^{2}-\frac{\alpha}{x_{2}^{*}n^{*}}(1-p_{S})(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{2}-x_{2}^{*})^{2}$
$\displaystyle+(r+\mu)[y_{0}+(1-p_{I})(y_{1}+y_{2})]\left\\{\frac{\alpha}{r+\mu}[\varepsilon_{0}+(1-p_{I})\varepsilon_{1}+(1-p_{S})(1-p_{I})\varepsilon_{2}]-1\right\\}$
$\displaystyle=$
$\displaystyle-\frac{\mu}{x_{0}^{*}}(x_{0}-x_{0}^{*})^{2}-\frac{\alpha}{x_{0}^{*}n^{*}}[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{0}-x_{0}^{*})^{2}$
$\displaystyle-\frac{\mu}{x_{1}^{*}}(1-p_{I})(x_{1}-x_{1}^{*})^{2}-\frac{\alpha}{x_{1}^{*}n^{*}}(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{1}-x_{1}^{*})^{2}$
$\displaystyle-\frac{\mu}{x_{2}^{*}}(1-p_{I})(x_{2}-x_{2}^{*})^{2}-\frac{\alpha}{x_{2}^{*}n^{*}}(1-p_{S})(1-p_{I})[y_{0}+(1-p_{I})(y_{1}+y_{2})](x_{2}-x_{2}^{*})^{2}$
$\displaystyle+(r+\mu)[y_{0}+(1-p_{I})(y_{1}+y_{2})](\mathcal{R}_{0}-1)$
$\displaystyle\leq$ $\displaystyle~{}0,$ (B.2)

which completes the proof.

##  Appendix C The existence and uniqueness of $\mathbf{\Psi}^{**}$ when
$\mathcal{R}_{0}>1$

Using the equations $\dot{y}_{0}=0$, $\dot{y}_{1}=0$, $\dot{y}_{2}=0$ in
$\dot{\mathbf{\Psi}}=\mathbf{0}$ to obtain $x_{i}^{**}$ as a function of
$y_{i}^{**}$, $i=0,1,2$. Then, substituting the results into the equations
$\dot{x}_{0}=0$, $\dot{x}_{1}=0$, $\dot{x}_{2}=0$,

$\left\\{\begin{aligned}
0=&~{}\varepsilon_{0}\Lambda-(r+\mu)y_{0}^{**}-\frac{\mu
n^{*}(r+\mu)y_{0}^{**}}{\alpha[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]},\\\
0=&~{}\varepsilon_{1}\Lambda-(r+\mu)y_{1}^{**}-\frac{\mu
n^{*}(r+\mu)y_{1}^{**}}{\alpha[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]},\\\
0=&~{}\varepsilon_{2}\Lambda-(r+\mu)y_{2}^{**}-\frac{\mu
n^{*}(r+\mu)y_{2}^{**}}{\alpha(1-p_{S})[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]}.\\\
\end{aligned}\right.$ (C.1)

In Eq. (C.1), we multiply the first equation by $\alpha/(r+\mu)$, the second
equation by $\alpha(1-p_{I})/(r+\mu)$, and the third equation by
$\alpha(1-p_{S})(1-p_{I})/(r+\mu)$:

$\left\\{\begin{aligned}
0=&~{}\frac{\alpha}{r+\mu}\varepsilon_{0}\Lambda-\alpha y_{0}^{**}-\frac{\mu
n^{*}y_{0}^{**}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})},\\\
0=&~{}\frac{\alpha}{r+\mu}(1-p_{I})\varepsilon_{1}\Lambda-\alpha(1-p_{I})y_{1}^{**}-\frac{\mu
n^{*}(1-p_{I})y_{1}^{**}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})},\\\
0=&~{}\frac{\alpha}{r+\mu}(1-p_{S})(1-p_{I})\varepsilon_{2}\Lambda-\alpha(1-p_{S})(1-p_{I})y_{2}^{**}-\frac{\mu
n^{*}(1-p_{I})y_{2}^{**}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})}.\\\
\end{aligned}\right.$ (C.2)

Summing up the three equations in Eq. (C.2) and using $n^{*}=\Lambda/\mu$ (see
Eq. (4)), we have

$\mathcal{R}_{0}-\alpha[y_{0}^{**}+(1-p_{I})y_{1}^{**}+(1-p_{S})(1-p_{I})y_{2}^{**}]-1=0.$
(C.3)

Therefore, to ensure
$y_{0}^{**}+(1-p_{I})y_{1}^{**}+(1-p_{S})(1-p_{I})y_{2}^{**}>0$, which is a
necessary condition for $y_{0}^{**}>0$, $y_{1}^{**}>0$, $y_{2}^{**}>0$, we
have $\mathcal{R}_{0}>1$. However, we have not yet proved that
$\mathcal{R}_{0}>1$ is a sufficient condition for $y_{0}^{**}>0$,
$y_{1}^{**}>0$, $y_{2}^{**}>0$.

According to Eq. (C.2), we can ensure $y_{0}^{**}>0$, $y_{1}^{**}>0$,
$y_{2}^{**}>0$ if we can confirm
$y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})>0$. We will try to illustrate the
opposite case, $y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})<0$, cannot happen.
Let us further write Eq. (C.2) as

$\left\\{\begin{aligned}
y_{0}^{**}=&~{}\dfrac{\dfrac{\alpha}{r+\mu}\varepsilon_{0}\Lambda}{\alpha+\dfrac{\mu
n^{*}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})}},\\\
y_{1}^{**}=&~{}\dfrac{\dfrac{\alpha}{r+\mu}\varepsilon_{1}\Lambda}{\alpha+\dfrac{\mu
n^{*}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})}},\\\
y_{2}^{**}=&~{}\dfrac{\dfrac{\alpha}{r+\mu}(1-p_{S})\varepsilon_{2}\Lambda}{\alpha(1-p_{S})+\dfrac{\mu
n^{*}}{y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})}}.\\\ \end{aligned}\right.$
(C.4)

Then, it can be judged that $y_{0}^{**}$ and $y_{1}^{**}$ have the same sign,
because the denominators are equal and the numerators are both positive. The
case $y_{0}^{**}<0$, $y_{1}^{**}<0$ is possible only if their denominators
$\alpha+\mu n^{*}/[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]<0$. In this
case, we have $\alpha(1-p_{S})+\mu
n^{*}/[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]<\alpha+\mu
n^{*}/[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]<0$, which means
$y_{2}^{**}<0$ as well. Then, we have
$y_{0}^{**}+(1-p_{I})y_{1}^{**}+(1-p_{S})(1-p_{I})y_{2}^{**}<0$ because
$y_{0}^{**}<0$, $y_{1}^{**}<0$, $y_{2}^{**}<0$, which is inconsistent with our
previous conclusion. Therefore, $y_{0}^{**}>0$, $y_{1}^{**}>0$ must hold.

The remaining question is the sign of $y_{2}^{**}$. According to the second
equation in Eq. (1), we have

$y_{0}^{**}=\frac{\alpha}{r+\mu}x_{0}^{**}[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]/n^{*}.$
(C.5)

Since we have $y_{0}^{**}>0$, we know
$x_{0}^{**}[y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})]>0$ as well. The sign
of $x_{0}^{**}$ can be easily judged: if $x_{0}=0$, then
$\dot{x}_{0}=\varepsilon_{0}\Lambda>0$ so that $x_{0}<0$ never happens if the
system starts evolving from a meaningful initial state where $x_{0}>0$.
Therefore, $x_{0}^{**}>0$ and $y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})>0$
must hold as well. Then, $y_{2}^{**}>0$ is ensured by the third equation in
Eq. (C.4).

Therefore, $\mathcal{R}_{0}>1$ is a sufficient and necessary condition for
$y_{0}^{**}>0$, $y_{1}^{**}>0$, $y_{2}^{**}>0$.

To check if the solution of $y_{0}^{**}$, $y_{1}^{**}$, and $y_{2}^{**}$
really exists, we can study the existence of
$y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})$. Then, the solution of
$y_{0}^{**}$, $y_{1}^{**}$, and $y_{2}^{**}$ can be naturally obtained by Eq.
(C.4). For convenience, we denote
$Y=y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})$. Multiplying the three
equations of $y_{0}^{**}$, $y_{1}^{**}$, $y_{2}^{**}$ in Eq. (C.4) by $1$,
$1-p_{I}$, $1-p_{I}$, and adding them together, we get

$Y=\dfrac{\dfrac{\alpha}{r+\mu}\varepsilon_{0}\Lambda}{\alpha+\dfrac{\mu
n^{*}}{Y}}+(1-p_{I})\dfrac{\dfrac{\alpha}{r+\mu}\varepsilon_{1}\Lambda}{\alpha+\dfrac{\mu
n^{*}}{Y}}+(1-p_{I})\dfrac{\dfrac{\alpha}{r+\mu}(1-p_{S})\varepsilon_{2}\Lambda}{\alpha(1-p_{S})+\dfrac{\mu
n^{*}}{Y}},$ (C.6)

which can be simplified as follows when $Y\neq 0$.

$aY^{2}+bY+c=0,$ (C.7)

where

$\left\\{\begin{aligned} a=&~{}\dfrac{\alpha}{\mu n^{*}}(1-p_{S}),\\\
b=&~{}1-p_{S}+1-\mathcal{R}_{0}+\frac{\alpha}{r+\mu}p_{S}(\varepsilon_{0}+(1-p_{I})\varepsilon_{1}),\\\
c=&~{}\frac{\mu n^{*}}{\alpha}(1-\mathcal{R}_{0}).\\\ \end{aligned}\right.$

We can see that $a>0$ always holds. The signs of $b$ and $c$, however, depend
on $\mathcal{R}_{0}$. Then, the simple use of Vieta’s theorem can judge the
existence of $Y$. When $\mathcal{R}_{0}<1$, we have $c/a>0$ and $-b/a<0$;
therefore, both roots are negative. When $\mathcal{R}_{0}>1$, we have $c/a<0$;
therefore, one root is positive and the other is negative. The positive root
is the unique solution of $Y$ and
$y_{0}^{**}+(1-p_{I})(y_{1}^{**}+y_{2}^{**})$. Then, Eq. (C.4) can provide the
unique solution of $y_{0}^{**}$, $y_{1}^{**}$, and $y_{2}^{**}$ (note that the
negative root of $Y$ cannot lead to positive $y_{0}^{**}$, $y_{1}^{**}$, and
$y_{2}^{**}$ because we have previously shown $y_{0}^{**}>0$, $y_{1}^{**}>0$,
and $y_{2}^{**}>0$ if $\mathcal{R}_{0}>1$).

## References

  * [1] Lynne Peeples.  “Face masks for COVID pass their largest test yet”.  Nature (Lond.), 2021.
  * [2] William Ogilvy Kermack and Anderson G McKendrick.  “A contribution to the mathematical theory of epidemics”.  Proceedings of the Royal Society of London A, vol. 115, no. 772, 700–721, 1927.
  * [3] Nathaniel S Barlow and Steven J Weinstein.  “Accurate closed-form solution of the SIR epidemic model”.  Physica D: Nonlinear Phenomena, vol. 408, 132540, 2020.
  * [4] Steven J Weinstein, Morgan S Holland, Kelly E Rogers and Nathaniel S Barlow.  “Analytic solution of the SEIR epidemic model via asymptotic approximant”.  Physica D: nonlinear phenomena, vol. 411, 132633, 2020.
  * [5] Chaoqian Wang and Chaochao Huang.  “An epidemic model with the closed management in Chinese universities for COVID-19 prevention”.  In Journal of Physics: Conference Series, volume 1707, page 012027\. IOP Publishing, 2020.
  * [6] Xinwei Wang, Haijun Peng, Boyang Shi et al.  “Optimal vaccination strategy of a constrained time-varying SEIR epidemic model”.  Communications in Nonlinear Science and Numerical Simulation, vol. 67, 37–48, 2019.
  * [7] Feng Fu, Daniel I Rosenbloom, Long Wang and Martin A Nowak.  “Imitation dynamics of vaccination behaviour on social networks”.  Proceedings of the Royal Society B: Biological Sciences, vol. 278, no. 1702, 42–49, 2011.
  * [8] Xinyu Wang, Danyang Jia, Shupeng Gao et al.  “Vaccination behavior by coupling the epidemic spreading with the human decision under the game theory”.  Applied Mathematics and Computation, vol. 380, 125232, 2020.
  * [9] Muntasir Alam, Kazuki Kuga and Jun Tanimoto.  “Three-strategy and four-strategy model of vaccination game introducing an intermediate protecting measure”.  Applied Mathematics and Computation, vol. 346, 408–422, 2019.
  * [10] Kazuki Kuga and Jun Tanimoto.  “Which is more effective for suppressing an infectious disease: imperfect vaccination or defense against contagion?”.  Journal of Statistical Mechanics: Theory and Experiment, vol. 2018, no. 2, 023407, 2018.
  * [11] Muntasir Alam, Masaki Tanaka and Jun Tanimoto.  “A game theoretic approach to discuss the positive secondary effect of vaccination scheme in an infinite and well-mixed population”.  Chaos, Solitons & Fractals, vol. 125, 201–213, 2019.
  * [12] Laijun Zhao, Qin Wang, Jingjing Cheng et al.  “Rumor spreading model with consideration of forgetting mechanism: A case of online blogging livejournal”.  Physica A: Statistical Mechanics and its Applications, vol. 390, no. 13, 2619–2625, 2011.
  * [13] Laijun Zhao, Wanlin Xie, H Oliver Gao et al.  “A rumor spreading model with variable forgetting rate”.  Physica A: Statistical Mechanics and its Applications, vol. 392, no. 23, 6146–6154, 2013.
  * [14] Chaoqian Wang.  “Dynamics of conflicting opinions considering rationality”.  Physica A: Statistical Mechanics and its Applications, vol. 560, 125160, 2020.
  * [15] Chaoqian Wang, Ziwei Wang and Qiuhui Pan.  “Injurious information propagation and its global stability considering activity and normalized recovering rate”.  Plos One, vol. 16, no. 10, e0258859, 2021.
  * [16] Cruz Vargas-De-León.  “On the global stability of SIS, SIR and SIRS epidemic models with standard incidence”.  Chaos, Solitons & Fractals, vol. 44, no. 12, 1106–1110, 2011.
  * [17] Jianquan Li, Yanni Xiao, Fengqin Zhang and Yali Yang.  “An algebraic approach to proving the global stability of a class of epidemic models”.  Nonlinear Analysis: Real World Applications, vol. 13, no. 5, 2006–2016, 2012.
  * [18] Sanusi Side, Wahidah Sanusi, Muhammad Kasim Aidid and Sahlan Sidjara.  “Global stability of SIR and SEIR model for tuberculosis disease transmission with lyapunov function method”.  Asian Journal of Applied Sciences, vol. 9, no. 3, 87–96, 2016.
  * [19] Hongbin Guo, Michael Y Li and Zhisheng Shuai.  “Global stability of the endemic equilibrium of multigroup SIR epidemic models”.  Canadian Applied Mathematics Quarterly, vol. 14, no. 3, 259–284, 2006.
  * [20] Ruoyan Sun.  “Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence”.  Computers & Mathematics with Applications, vol. 60, no. 8, 2286–2291, 2010.
  * [21] Yoshiaki Muroya, Yoichi Enatsu and Toshikazu Kuniya.  “Global stability for a multi-group SIRS epidemic model with varying population sizes”.  Nonlinear Analysis: Real World Applications, vol. 14, no. 3, 1693–1704, 2013.
  * [22] Lefeng Cheng, Linfei Yin, Jianhui Wang et al.  “Behavioral decision-making in power demand-side response management: A multi-population evolutionary game dynamics perspective”.  International Journal of Electrical Power & Energy Systems, vol. 129, 106743, 2021.
  * [23] Lefeng Cheng, Yang Chen and Guiyun Liu.  “2PnS-EG: A general two-population $n$-strategy evolutionary game for strategic long-term bidding in a deregulated market under different market clearing mechanisms”.  International Journal of Electrical Power & Energy Systems, vol. 142, 108182, 2022.
  * [24] Weiqiang Li, Jin Zhou and Jun-an Lu.  “The effect of behavior of wearing masks on epidemic dynamics”.  Nonlinear Dynamics, vol. 101, no. 3, 1995–2001, 2020.
  * [25] João AM Gondim.  “Preventing epidemics by wearing masks: An application to COVID-19”.  Chaos, Solitons & Fractals, vol. 143, 110599, 2021.
  * [26] Pierre Auger and Ali Moussaoui.  “On the threshold of release of confinement in an epidemic SEIR model taking into account the protective effect of mask”.  Bulletin of Mathematical Biology, vol. 83, no. 4, 1–18, 2021.
  * [27] Nurudeen O Lasisi and Kolawole A Adeyemo.  “Modeling the effect of distancing and wearing of face masks on transmission of Covid-19 infection dynamics”.  Journal of Complexity in Health Sciences, vol. 4, no. 1, 10–20, 2021.
  * [28] Lili Han, Qiuhui Pan, Baolin Kang and Mingfeng He.  “Effects of masks on the transmission of infectious diseases”.  Advances in Difference Equations, vol. 2021, no. 1, 1–17, 2021\.
  * [29] Scott E Page.  The model thinker: What you need to know to make data work for you.  Basic Books, 2018.
  * [30] Serge. Galam.  Sociophysics: A physicist’s modeling of psycho-political phenomena.  Springer, 2016.
  * [31] Pauline Van den Driessche and James Watmough.  “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission”.  Mathematical Biosciences, vol. 180, no. 1-2, 29–48, 2002.
  * [32] Joseph P La Salle.  The stability of dynamical systems.  SIAM, 1976.
  * [33] Mohammed A Boraey.  “An analytical model for the effective filtration efficiency of single and multiple face masks considering leakage”.  Chaos, Solitons & Fractals, vol. 152, 111466, 2021.