Understanding the Role of Environmental Transmission on COVID-19 Herd Immunity and Invasion Potential

COVID-19 is caused by the SARS-CoV-2 virus, which is mainly transmitted directly between humans. However, it is observed that this disease can also be transmitted through an indirect route via environmental fomites. The development of appropriate and effective vaccines has allowed us to target and anticipate herd immunity. Understanding of the transmission dynamics and the persistence of the virus on environmental fomites and their resistive role on indirect transmission of the virus is an important scientific and public health challenge because it is essential to consider all possible transmission routes and route specific transmission strength to accurately quantify the herd immunity threshold. In this paper, we present a mathematical model that considers both direct and indirect transmission modes. Our analysis focuses on establishing the disease invasion threshold, investigating its sensitivity to both transmission routes and isolate route-specific transmission rate. Using the tau-leap algorithm, we perform a stochastic model simulation to address the invasion potential of both transmission routes. Our analysis shows that direct transmission has a higher invasion potential than that of the indirect transmission. As a proof of this concept, we fitted our model with early epidemic data from several countries to uniquely estimate the reproduction numbers associated with direct and indirect transmission upon confirming the identifiability of the parameters. As the indirect transmission possess lower invasion potential than direct transmission, proper estimation and necessary steps toward mitigating it would help reduce vaccination requirement.


Introduction
Coronaviruses are enveloped RNA viruses that use mammals and birds as hosts and have the ability to cause various types of respiratory symptoms (Wardeh et al. 2021;Zhu et al. 2020;. Two distinguished strains of this virus, namely SARS-CoV and MERS-CoV, have caused several epidemic outbreaks during the last two decades at several places around the world (Zhu et al. 2020). The ubiquity of this virus along with its large genetic diversity and increasing animal-human interactions has amplified the likelihood of the emergence of a coronavirus infection (Huang and Wang 2021). The most recent outbreak of the virus was caused by the novel strain SARS-CoV-2 that led to the recent pandemic of the coronavirus disease 2019 .
In the last two years, significant improvement has been done in understanding the transmission routes and pathways of COVID-19 (Azuma et al. 2020;Rothe et al. 2020;Yu and Yang 2020;Morawska et al. 2020;Pitol and Julian 2021;Castaño et al. 2021). The onset of this disease is usually characterized by symptoms like fever, cough, and sore throat, and in some cases, the severity of the disease leads to shortness of breath. Virus particles discharged through nostrils and mouth during breathing, talking, sneezing, and/or coughing may transmit the disease to other host. COVID-19 patients may spread the disease at least 1-3 days before the onset of their symptoms (Wormser 2020). Furthermore, in many cases (17. 8% Mizumoto et al. 2020, 30.8% Nishiura et al. 2020, it has been shown that patients tend to be asymptomatic or simply develop very mild symptoms throughout the entire infectious period. Consequently, patients who are infectious and transmit the disease may go unnoticed, which can be a key driver that undermines any efforts to contain the disease (Bai et al. 2020;Rothe et al. 2020). Another potential driver for transmission could be the prolonged sustenance of the virus on environmental fomites (Vardoulakis et al. 2020;Azuma et al. 2020;Pitol and Julian 2021). In experimental setup, SARS-CoV-2 was found stable on plastic and stainless steel up to 72 h (van Doremalen et al. 2020). On plastic and human skin surfaces, variants of SARS-CoV-2 maintained infectivity for several hours (Hirose et al. 2020(Hirose et al. , 2022. Gidari et al. reported infectious existence of this virus on plastic and glass for more than 120 h and on stainless steel for more than 72 h (Gidari et al. 2021). Infectious virus was detected even after 7 days on a sample of surgical masks (Chin et al. 2020). SARS-CoV-2 survival for up to 1, 5, and 10 days was reported on fake fur, plastic, and mink fur, respectively (Brown et al. 2021). In artificial saliva, it was found stable for at least 90 min (Smither et al. 2020). Live SARS-CoV-2 RNA was detected on 8.3% of the high-touch surfaces in the public locations during a COVID-19 outbreak in Massachusetts (Harvey et al. 2021). The above literature suggests that an additional key driver for COVID-19 outbreak could be the prolonged sustenance of the virus on environmental fomites. However, the effectiveness of surface disinfection is highly dependent on the prevalence and the frequency of contact as well as environmental conditions (Gidari et al. 2021;Pottage et al. 2021;Wilson et al. 2021). For instance, approximately 30% of disease transmissions on the Diamond Princess cruise ship were reported to be related to fomite-mediated transmission (Azimi et al. 2021), whereas in China this percentage is reported to be 45-62% (Yang and Wang 2021). In hospital setting, 27% of the environmental surfaces were reported to contain SARS-CoV-2 RNA even though disinfectant were sprayed twice . Further details pertaining to the deposition, survival, and transmission of the virus can be accessed in Leung (2021), Castaño et al. (2021), Aydogdu et al. (2021), and Gonçalves et al. (2021). The environmental transmission has also been observed to play a critical role in the persistence and inter annual epidemics for other communicable diseases (Vergara-Castaneda et al. 2012;Lopman et al. 2012;McKinney et al. 2006;Breban et al. 2009;Al-Tawfiq and Memish 2016).
Environmental route of transmission has been modeled mathematically for several infectious diseases and was proved to hold important implications for disease control (Eisenberg et al. 2005;Zhao et al. 2012). For instance, environmental transmission modulates the periodicity in avian influenza outbreak Rohani et al. 2009;Wang et al. 2012). It is also associated with spatial diffusion of avian influenza (Li et al. 2019). Recently, several mathematical models have been proposed regarding fomite-mediated transmission of COVID-19 (Yang and Wang 2020;Stutt et al. 2020;Yang and Wang 2021;Wijaya et al. 2021;Rwezaura et al. 2021). However, the role of fomite-mediated transmission in crucial public health issues, such as herd immunity, has yet to be substantially explored. Moreover, it would be of interest to evaluate which transmission pathway has a higher invasion potential. In this study, we have classified the transmission routes into two types-direct and indirect. Direct transmission refers to transmission of the infection that comes directly from an infectious to a susceptible individual. In contrast, indirect transmission refers to the deposition of the virus particles by an infectious individual on environmental fomites followed by inoculation of the virus by a susceptible individual who, in turn, becomes infectious.
The remaining part of the paper is structured as follows. We initially present the mathematical model in Sect. 2. We then analyze our model to establish the disease invasion threshold and investigate its sensitivity to both transmission routes in Sect. 3.1. Section 3.2 presents the stochastic simulation wherein we investigate the invasion potential of both transmission routes. Consequently, we fit our model to the early epidemic data obtained from several countries, along with identifiability analysis to quantify route-specific transmission strength (in Sect. 3.3), thereby measuring the vaccination requirement according to initial outbreak data for the acquisition of herd immunity, which is presented in Sect. 3.4. Finally, we discuss and summarize our findings in Sect. 4.

Epidemic Model
We divided the human host into four different compartments: susceptible (S), infectious ( A), confirmed infected (I ), and recovered (R). Furthermore, F represents viruses on environmental carriers or fomites (Fig. 1). Susceptible individuals are the ones who can contract the disease following exposure to the virus. Once a susceptible person is exposed to the virus through direct or indirect contact with an infectious agent, they may either become infectious and express symptoms after a latency period or may not express symptoms albeit transmit the disease. Depending on symptom expression, public health guidelines and the capacity of public health authority to isolate infectious individuals, the infectious individuals may be confirmed/identified as infected, or may remain unnoticed and remain infectious. For simplicity, we assume the confirmed infected individuals no longer transmit the disease. An infected individual may remain infectious throughout his/her whole infectious lifetime and pass through S − → A − → R pathway; or an infected individual may remain infectious in first few days until s/he becomes confirmed at some point of his/her infectious lifetime and pass through S − → A − → I − → R pathway.
When infectious individuals talk loudly, cough, or sneeze, numerous virus particles exit from their respiratory organs and can be deposited on surfaces in the environment, where they can survive for a long time (van Doremalen et al. 2020;Hirose et al. 2020Hirose et al. , 2022Gidari et al. 2021;Chin et al. 2020;Brown et al. 2021;Smither et al. 2020;Harvey et al. 2021) and be carried away by a new host afterward. Virus particles deposited on environmental fomites belong to the F compartment. The rate at which asymptomatic individuals deposit the virus on fomites is and the viral particles decay naturally at a rate of ξ . We assume that the viral population is large enough to describe the viral population dynamics by the following two processes: deposit and decay. In addition, we consider that the probability of infection from a virus picked up from the environmental fomites is a function of daily viral exposure to the environmental fomites. Precisely, we assume that a constant fraction ρ of virions (F) is picked up by each susceptible individual per day and may cause infection with a probability of g(F d ), where F d = ρ F, which reflects the daily pick up rate. We consider the following two types of functional forms for g, which are written as a function of F only for simplicity, as ρ is assumed constant.
where α = πρ. Under these assumptions, we obtain the following system of differential equations

Invasion Threshold
The invasion threshold of the disease is determined by the existence and stability of the equilibria.

Proposition 1
The model has a unique disease-free equilibrium (DFE), E 0 . In addition, it has an endemic equilibrium (EE), E which exists for R 0 > 1.
The basic reproduction number, R 0 , is a crucial threshold for characterizing the dynamics of an outbreak. It refers to the average number of secondary infections caused by the introduction of one infectious individual in a completely susceptible population. Here, we designate A and F as the diseased class. The DFE is given by E 0 = μ , 0, 0, 0, 0 . The next generation matrix at E 0 is given by (Please refer to "Appendix A" for details) where k A = ω + γ A + μ and k I = γ I + μ + δ. K k, j provides the expected number of secondary infections in class k produced by a single incident in class j. Recent studies (Yang and Wang 2020;Stutt et al. 2020;Yang and Wang 2021;Wijaya et al. 2021), K 2,1 is considered 0 without considering the environmental fomites (F) as infectious. In this study, we consider both A and F as infectious compartments, and we interpret K 1,1 as the number of new secondary infectious individuals caused by one infectious individual during his/her entire infectious period, K 2,1 as the number of virus particles spread by an infectious individual throughout his/her entire infectious period, and K 1,2 as the number of new infected individual caused by one virus particle in the environment throughout its entire active period. The spectral radius of K is the basic reproduction number R 0 , which is given below The basic reproduction number can be rearranged and expressed in the following form: Here, we can highlight the role of direct and indirect transmission. At E 0 , an asymptomatic infectious individual transmits disease to μ β individuals per day and the total infectious period is 1 k A days. Therefore, R 0H = μ β 1 k A is the expected number of new infections generated from an infected individual throughout his/her entire infectious period. In contrast, infectious individuals deposit virus particles on environmental fomites at a rate of . The total virus particles deposited in the environment throughout the entire infectious period of a single infected individual is k A . Each virus particle can subsequently infect a person with a probability of α in Case I and 1 − e −α ≈ α in Case II per day. A virus particle can survive on environmental fomites for an expected duration of 1 ξ days. Hence, the expected number of infected individuals caused by a single virus particle in a completely susceptible environment is μ 1 ξ α. Therefore, R 0F = πρ μξ k A quantifies the expected number of secondary infections as a result of indirect transmission.
A comparison of the role of direct and indirect transmission is portrayed in Fig. 2, which shows that R 0 increases linearly with R 0H . In contrast, R 0 increases faster than linear with R 0F until R 0 < 1. When R 0 > 1, the impact of indirect transmission diminishes as R 0 increases, i.e., indirect transmission plays a crucial role if R 0 is near 1.

Theorem 1
The DFE is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1.
Please refer to "Appendix B" for detailed proof. The EE is given by and for Case II, A * is given by 1], otherwise the system will be unbounded below. Hence, we obtain which has, at most, one positive solution for A * as 1−(R 0H +1)μ β A * ≥ 0 since μ << 1. For Case I, the expression of A * demonstrates that the EE exists for R 0H + R 0F ≥ 1. However, it is cumbersome to deduce the condition for the existence of A * in Case II. Therefore, we chose to use the center manifold theorem (Theorem 4.1 in Castillo-Chavez and Song 2004), which could help us characterize the existence and nature of the EE near R 0 = 1 using the Jacobian at the DFE for both the cases simultaneously.

Theorem 2
The model exhibits forward bifurcation at R 0 = 1.
Please refer to "Appendix C" for detailed proof. The above analysis nominates R 0 = 1 as the disease invasion threshold. The expression of R 0 clarifies the role of both direct and indirect transmissions. R 0 can be used to measure the control efforts required to mitigate or stop the spread of the disease. However, if the infection is carried by more than one types of host, the use of R 0 leads to a distinct underestimation of the requirements (Bani-Yaghoub et al. 2012;Pauline 2017). In such cases, the type reproduction number provides a significantly more accurate estimation of the required control efforts (Roberts and Heesterbeek 2003;Heesterbeek and Roberts 2007).

Type Reproduction Number
It is essential to gain a clear understanding of the explicit role of human carriers (direct transmission) and environmental carriers (indirect transmission) in spreading the virus so as to decide feasible and effective control strategies. The basic reproduction number properly defines the invasion threshold, but this number cannot distinguish the pathway-specific transmission. In this section, we use the concept of type reproduction number (Roberts and Heesterbeek 2003;Heesterbeek and Roberts 2007) to investigate the pathway-specific transmission. Type reproduction number for host, T H , refers to the expected number of infectious individuals caused by one infectious individual in a completely susceptible environment, either by direct or indirect transmission. Following the notation in Roberts and Heesterbeek (2003), let I 5 be the 5 × 5 identity matrix, P H = [ph i j ] be the projection matrix defined by, ph 11 = 1, and ph i j = 0 when i = 1 or j = 1. E H is the unit column matrix with its first element equal to 1. Then, and we have the following properties (Heesterbeek and Roberts 2007).
This reduction can be achieved by means of vaccinating susceptible individuals (as this will reduce the number of available susceptible individuals) or by quarantining infectious individuals (as this will reduce their infectious period).
Therefore, the invasion threshold can be refined in terms of the type reproduction number as T H = 1. Furthermore, although the expression for R 0 is difficult to interpret from the biological point of view, the expression for T H can be easily understood as the total number of the expected secondary infected individuals as a result of both the direct and indirect transmissions caused by one infectious individual in a completely susceptible environment. Figure 3 shows a comparison between R 0 and T H illustrating that they both coincide at 1, but R 0 > T H below 1 and R 0 < T H above 1, which may also be confirmed by using simple algebra. Rephrasing the invasion threshold not only allows us to provide a biological interpretation but also leads us to differentiate the pathway-specific transmission strength and infer the relative requirement of the subsequent control measures (Heesterbeek and Roberts 2007). If the value of T H is known, we can estimate the vaccination requirement. Further, if we can distinguish R 0H and R 0F , i.e., isolate pathway-specific transmission strength, we will be able to estimate required strictness in maintaining quarantine measures, and requirement of cleanliness and maintaining personal hygiene to an appropriate degree.

Stochasticity in Invasion
The expression of invasion threshold, T H demonstrates that it is equally sensitive to R 0H and R 0F . For T H slightly greater than 1, there might exist a nonzero probability of disease extinction. Direct transmission depends on one successful transmission from one host to another, whereas indirect transmission hinges on two successful transmissions-one from the original host to fomite and then back from fomite to another host upon survival. To inspect the potential impact of stochasticity associated with these different transmission pathways on the invasion potential, we performed a stochastic simulation using the Modified Poisson Tau-Leap algorithm (Cao et al. 2005). The technique has been explained in "Appendix D" and parameter values have also been presented. To understand the invasion potential, we simulated our model for a duration of 1 year for T H = R 0H + R 0F = 1.1, 1.2, where both R 0H and R 0F vary from 0% to 100% of T H to maintain the specified value of T H . We ran 1000 simulations for each case. Among the 1000 simulations, the fraction of number of times infectious individuals, A(t) that reaches zero provides us an approximate extinction probability, which is plotted in Fig. 4. The figure shows that when T H = 1.1, R 0H ≤ 0.88 and R 0F ≥ 0.22, the disease goes extinct by the end of 1 year. In contrast, when R 0H > 0.88 and R 0F < 0.22, the extinction probability decreases to approximately 0.9 for both cases I (g 1 (F)) & case II (g 2 (F)). In contrast, when T H = 1.2, R 0H ≤ 0.84 and R 0F ≥ 0.36, the disease goes extinct by the end of 1 year. However, when R 0H > 0.84 and R 0F < 0.36, the extinction probability decreases to approximately 0.86 for both cases I and II. Therefore, the extinction probability decreases with increasing T H and indirect transmission has a higher chance of extinction compared to direct transmission.  Fig. 4 Extinction probability. Approximate extinction probabilities for different combinations of R 0H and R 0F corresponding to the same T H have been shown. The figure on the left is for 0 ≥ R 0F , R 0H ≥ 1.1, whereas the figure on the right is for 0 ≥ R 0F , R 0H ≥ 1.2 (color figure online)

Vaccination Threshold for Herd Immunity
In transmissible diseases, viral pathogens in existing hosts attempt to find another host to survive, proliferate and complete the cycle of transmission. At this stage, should the susceptible hosts become significantly scarce, which means that the virus is unable to find a suitable host for transmission, the transmission cycle breaks and the virus goes extinct. This is possible if the overall population has a sufficient number of immune individuals, which is defined as the state of herd immunity (Rasmussen 2020;McDermott 2021). It is a dynamic threshold that depends on the reproduction number and, consequently, on the disease transmission rate. In our present problem, this threshold is v c = 1 − 1 T H , i.e., if v c fraction of host becomes immune to the virus by vaccination or recovering from the infection, the pandemic will end. From the expression of v c , it is comprehensible that if the reproduction number increases (for instance, as a consequence of increasing the transmission rate), the herd immunity threshold increases as well. Therefore, it is not rational to define a rigid threshold v c that is less than unity, and thus remove all preventive measures.
One important information that T H provides us with is the transmission pathway specific requirement of control measures for herd immunity. In the COVID-19 case, this allows us to distinguish between the requirements of control measures against direct transmission and those against indirect transmission. We can minimize indirect transmission in T H by conforming to safety practices, such as general cleanliness, good hygiene, and disinfecting surfaces, while vaccination and different forms of quarantine measures can reduce the direct transmission in T H . If we do not use any measures to prevent environmental transmission, the vaccination requirement for herd immunity, v c , would be v c,max = 1 − 1 R 0H +R 0F . Further, if we take measures for reducing the environmental transmission, the threshold (v c ) would then satisfy the inequality 1− 1 R 0H +R 0F > v c > 1− 1 R 0H . Provided that the environmental transmission could be completely stopped, the vaccination threshold would be reduced to v c,min = 1 − 1 R 0H .
However, it is challenging to isolate the correct pathway specific transmission strength as fitting a model with non-identifiable set of parameters may induce errors of attributing the contribution of one transmission pathway to another. Therefore, we first check the identifiability of the parameters, which confirms the uniqueness of the estimated parameter values, and thus isolates the pathway specific transmissibility. To clarify this with examples, we fit our model with daily active cases of early epidemic data from Nigeria, Bangladesh, and USA (considering g 1 (F) as the fomite to human transmission function), and we then estimate the parameters β and α. The daily active cases data are taken from worldometer (https://www.worldometers.info/ coronavirus/).

Identifiability and Fitting
Let us assume, X = (S, A, I , R, F) and denote the right side of the system (1) by F. Further, P = (β, α) be the vector of parameters to be estimated. We assume the remaining parameters to be known and summarize in Table 1 with proper citation. Here, I (t, P) is the vector of observable and i(t, P) is the observed data at t = 1, 2, . . . , 40 days. We assume i(t, P) follow Poisson distributed with mean I (t, P), then the maximum likelihood function will be: As ln is a monotonically increasing function, we minimize the negative log likelihood function (NLF) instead of maximizing the likelihood function for computational convenience. The NLF is reduced to: As the last term in the above equation remains unchanged, it is sufficient to minimize the sum of the first two terms. Therefore, the fitting process reduces to a minimization problem as, subject to d dt X (t, P) = F(X , P, t)  (2020) The above fitting problem will provide practically feasible and unique parameters values if P is identifiable. The parameters P is structurally identifiable if a unique solution X (t, P) exists for each P and a fixed initial condition. First, we estimate the fisher information matrix (F I M) and then compute the profile likelihoods to confirm the identifiability of the parameters.
We have observations at 40 distinct times, a system of 5-state variables, and two unknown parameters. Therefore, the sensitivity matrix M consists of 5 time-dependent 5 × 2 blocks A(t k ) . . .
The 2 × 2 FIM is F I M = M T M, which has 2 columns. Let us denote the parameter estimates asβ andα. We approximate the FIM numerically by perturbingβ to the valuesβ + = (1 + 0.001)β andβ − = (1 − 0.001)β, for which we integrate the model for each observation time. Then, we approximate the derivatives, fixed. This provides the first column. We repeat the same process forα to obtain the second column. Then, we check the rank of the matrix F I M, which is 2, which ensures that the parameters have no implicit dependency. This confirms the structural identifiability numerically. Further. we investigate practical identifiability to confirm whether the parameters estimated by fitting this model with this set of data are capable of differentiating the role of the different transmission pathways.
To investigate practical identifiability, we compute the profile likelihood of the parameters β and α. Profile likelihood reveals the dependency of the N L F on each parameter, and exposes the minimization of the N L F at the estimated value. The desired profile likelihoods are as follows: where β ∈ [β(1 − 0.05),β(1 + 0.05)] and α ∈ [α(1 − 0.05),α(1 + 0.05)]. Figure 5 shows the fitting along with the profile likelihood of the parameters which shows unique minima of the NLF at the estimated value of the parameters (second and third column) and hence confirms the identifiability which informs pathway specific transmission potential. The corresponding estimates of the reproduction numbers along with bounds for vaccination threshold are shown in Table 2 Figure 6 depicts the vaccination thresholds for these three different countries as a function of R 0F , and it clearly shows that the vaccination threshold would be decreased to a minimum value v c,min = 1 − 1 R 0H when R 0F = 0, i.e., the vaccination requirement would reach its minimum value v c,min , which is the y-intercept, if we manage to take sufficient measures to ensure no environmental/indirect transmission. In contrast, when the environmental transmission is partially halted, or if no measures are taken whatsoever, the vaccination threshold would be increased to a maximum value of v c,max = 1 − 1 R 0H +R 0F depending on the size of R 0F . Moreover, as the indirect transmission possess higher chance of extinction than direct transmission, the additional indirect transmission would not increase the probability of invasion. According to our estimation from early epidemic data, provided that we limit indirect transmission, the vaccination threshold can be reduced to a minimum of 0.191, 0.130, and 0.353, in the cases of Nigeria, Bangladesh, and USA, respectively. Table 1 summarizes other parameter values used in the simulations for these three countries. It is noteworthy that as time will pass by, the estimate of both T H and v c will change.

Discussion and Conclusion
SARS-CoV-2 can survive on different types of surfaces and has the potential to be transmitted to susceptible individuals. Therefore, our model has considered two types of transmission routes: human-to-human (direct transmission) and human to environmental fomites and then back to human (indirect transmission). Both transmission routes contribute to the reproduction number, and the degree of this epidemic is enhanced when the sum of the contribution of both direct and indirect transmissions in the type reproduction number exceeds one. The deterministic result shows that the invasion threshold is equally sensitive to both transmission routes. However, the stochastic simulation reveals that it is the indirect transmission that has a lower invasion potential compared to direct transmission.
Our analysis demonstrates that to develop effective control strategies, it is important to differentiate the role of these two different routes. The explicit modeling of both transmission routes and the estimation of the associated reproduction rates can allow us to gain a greater understanding of the indirect transmission epidemic potential and extent of efforts and measures that should be implemented in terms of disinfecting our proximal environment and maintaining personal hygiene. The analysis shows that the epidemic may persist even if direct transmission is reduced to 0 (for example, by social distancing and/or vaccination), whereas the reproduction number due to indirect transmission is > 1 (for example, due to lack of personal hygiene). Similarly, the epidemic may persist even if the indirect transmission is reduced to 0, whereas the reproduction number due to the direct transmission is > 1. It should be noted that it might not be practically feasible to reduce transmission from either routes to 0. If the reproduction number due to the direct transmission is < 1 but the type reproduction number is > 1, the epidemic could simply be terminated by maintaining strict cleanliness only. Moreover, the environmental transmission has lower invasion potential then the direct transmission. This nourishes the conclusion that, once the strength of the environmental transmission is known, direct transmission can be contained by using focused controls, such as the vaccination and/or different forms of quarantine.
Having obtained the sensitive quantification of the epidemic potential of the environmentally mediated transmission, mitigating environmental transmission by cleanliness, personal hygiene, and disinfection of the contaminated surfaces would reduce the requirement of human oriented control efforts. Note that individuals with either vaccine or acquired immunity may not be responsible for direct transmission but may play a plausible role in indirect transmission by acting as carrier. Besides, the herd immunity threshold is not a steady state; instead, it may increase due to the increasing transmission rate or increase in susceptible individuals due to loss of immunity. Therefore, having a fraction of individuals with immunity does not allow us to abort all preventive measures. Moreover, besides the implementation of vaccines, SARS-CoV-2 is rapidly developing mutations and it is likely that the vaccine-related antibodies may become ineffective to the new strains. We, therefore, conclude that all possible transmission routes need to be carefully considered and measured while vaccinating the population until the transmission is fully under control or declared eradicated.
Apart from fomite-mediated transmission, indirect transmission also includes transmission through aerosol particles deposited in the air by droplets, as suggested by Leung (2021), Castaño et al. (2021) and Aydogdu et al. (2021). However, aerosol mediated transmission is highly dependent on the respective air circulation and the physical structure of the venue. The deposition and decay rates can vary significantly. Furthermore, consideration of both indirect routes would lead to difficulties in terms of estimating the parameters associated with them. Despite this limitation, our study provides a clear indication that all possible transmission routes need to be carefully considered, and their transmission potential needs to be accurately quantified to mea-sure the required accurate threshold for achieving herd immunity. Lastly, in estimating herd immunity and planning control strategies, immunity loss and the probability of reinfection should be considered as well, two issues that form the core of our future research studies.
Funding This work is supported by the National Research Foundation of Korea (NRF) Grant funded by the Korea Government (NRF-2017R1E1A1A03069992).

Declarations
Conflict of Interest All the authors declare that there is no conflict of interest.
Ethical Approval Not applicable.

Consent for Publication
All the authors give consent.

Code Availability Available on request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Next Generation Matrix
To find the basic reproduction number, we follow the next generation method (van den Driessche and Watmough 2002) and system (1) is written in the following form: where F, V − , and V + are as follows: Here, the rate of appearance of new infections is represented by the elements of F and the rate of transfer of individuals into and out of the compartments are in V + and V − , respectively. The Jacobian of F at E 0 takes the following form.
As only the first 2 × 2 block of the above matrix is nonzero, whatever multiplied with it would result in a similar matrix with only a first 2 × 2 nonzero block. So, it suffices to consider a 2 × 2 matrix instead of a 5 × 5 matrix. Following van den Driessche and Watmough (2002), we write the sub-matrices F and V as

Appendix B: Proof to Theorem 1
Proof Conferring the symbols from equation (A1), we verify the following at E 0 . • Therefore, theorem 2 from van den Driessche and Watmough (2002) completes the proof.

Appendix C: Proof to Theorem 2
Proof For computational convenience, we assume k A as the bifurcation parameter. At R 0 = 1, The largest eigenvalue of the Jacobian of the system at DFE and k * A is 0 with multiplicity one. In addition, there is a nonnegative left eigenvector and a nonnegative right eigenvector corresponding to the 0 eigenvalue. The local dynamics are then determined by the two constants a and b defined in Theorem 4.1 of Castillo- Chavez and Song (2004). For our model, a = − β 2 and b = ξ 2 μ πρ k * 2 A , which indicates that a positive EE exists and it is locally asymptotically stable for R 0 > 1.
Here, (.) 2 stands for the element-wise square and * stands for the matrix product. Therefore, following Cao et al. (2005)  3: Initialize the current time t = 0; 4: At time t for state x, evaluate the propensity matrix a, its gradient ∇a, and a 0 = 12 j=1 a j ; 5: Find the set of indices, c i of the currently critical reactions for which a j (x(t)) > 0 and L j ≤ n c ; 6: Compute the largest possible time step τ that would not allow any propensity function to change its value by more than εa 0 , where the index j runs over the currently critical reactions identified in step 5. If there are no critical reactions assign τ = ∞. 7: If τ < 10 a 0 , reject it and run SSA (Gillespie 1976) and proceed to the next time step. Otherwise, continue to the next step. 8: Compute the sum of the propensity functions of the critical functions a c 0 (x(t)) = j∈c i a j (x) and generate τ = 1 a c 0 (x(t)) ln( 1 r ) where r ∼ U (0, 1). 9: If τ < τ , assign τ ← τ , for all the critical reactions R j , set k j = 0. 10: If τ ≥ τ , assign τ ← τ , generate j c as a sample of the integer random variable following j c = min j c ∈c i { j c | j c k=1 a k (x) > r 2 k∈c i a k (x)}. Set k j c = 1, and for all the other critical reactions set k j = 0, where j ∈ c i . 11: For all the noncritical reactions R j , generate k j ∼ Poisson(a j (x(t))τ ).
12: If min{x + 12 j=1 k j v j } < 0 assign τ ← τ 2 and return to step 9. 13: Assign t ← t + τ and x ← x + 12 j=1 k j v j . 14: Record (t, x) and return to step 5 until t reach the end of the time span. Furthermore, L j = min (v i j <0) i∈{1,2,...,5} x i |v i j | . Here, the square bracket represents the "greatest integer" operator. The simulation is completed according to Algorithm 1.