Is the Increased Transmissibility of SARS-CoV-2 Variants Driven by within or Outside-Host Processes?

Abstract

Understanding the factors that increase the transmissibility of the recently emerging variants of SARS-CoV-2 can aid in mitigating the COVID-19 pandemic. Enhanced transmissibility could result from genetic variations that improve how the virus operates within the host or its environmental survival. Variants with enhanced within-host behavior are either more contagious (leading infected individuals to shed more virus copies) or more infective (requiring fewer virus copies to infect). Variants with improved outside-host processes exhibit higher stability on surfaces and in the air. While previous studies focus on a specific attribute, we investigated the contribution of both within-host and outside-host processes to the overall transmission between two individuals. We used a hybrid deterministic-continuous and stochastic-jump mathematical model. The model accounts for two distinct dynamic regimes: fast-discrete actions of the individuals and slow-continuous environmental virus degradation processes. This model produces a detailed description of the transmission mechanisms, in contrast to most-viral transmission models that deal with large populations and are thus compelled to provide an overly simplified description of person-to-person transmission. We based our analysis on the available data of the Alpha, Epsilon, Delta, and Omicron variants on the household secondary attack rate (hSAR). The increased hSAR associated with the recent SARS-CoV-2 variants can only be attributed to within-host processes. Specifically, the Delta variant is more contagious, while the Alpha, Epsilon, and Omicron variants are more infective. The model also predicts that genetic variations have a minimal effect on the serial interval distribution, the distribution of the period between the symptoms’ onset in an infector–infectee pair.

1. Introduction

Since December 2020, the genomic surveillance effort in many countries has led to the detection of numerous variants of SARS-CoV-2 [1,2,3]. Some of these variants have exhibited increased transmissibility and raised concerns among public health authorities worldwide due to the risk that they will spread faster than vaccine production and distribution [2,4,5]. Understanding the mechanism that enhances the transmissibility of these variants is an essential step in devising methods to control their transmission [6].

The transmission of a virus between individuals is a complex multiscale process and consists of both within-host, and outside-host processes [3,7]. Many studies have provided information regarding the increased transmissibility of the dominant SARS-CoV-2 variants. Some of these have dealt with virus attributes that are associated with within-host processes, such as the viral load [8,9], the ability to initiate infection with a low dose of inocula [10], and the affinity of the virus spike receptor binding domain (RBD) to the host’s angiotensin-converting enzyme 2 (ACE2) [11,12]. Other studies have examined virus characteristics that influence outside-host processes, such as the stability of virus variants in aerosols [13,14] and their persistence on surfaces [15]. Despite these studies’ important information, they all suffer from two main limitations. First, they deal with a small, specific set of variants, while more transmissible variants continue to emerge. Second, as these studies examine a single attribute, it is impossible to compare these attributes’ contribution to the phenomenon of virus transmission as a whole.

Computational models can overcome these shortcomings, as they can cope with limited data, produce synthetic scenarios, and quantify the relative contribution of different parameters to the virus transmission [16,17,18,18]. While various models were used to examine the impact of highly transmissible variants, e.g., [5,19,20,21], very few models have tried to investigate possible reasons for the transmissibility of an emerging variant [20,21]. As these models describe the viral transmission among large populations, at the scale of a region or a major city, they lack details of transmission mechanisms between two individuals, especially those that occur in indoor spaces [22]. Therefore, such models can provide only a partial explanation for the influence of variants on the different aspects of the transmission process.

This work aims to determine whether the increased transmissibility of the SARS-CoV-2 variants (such as the Alpha, Epsilon, Delta, and the recently emerging Omicron) is the result of an enhanced within-host or outside-host behavior. To answer that question, we utilized a validated, mechanistic, multiscale, mathematical model to describe the transmission process of SARS-CoV-2 in a household. This model describes the specific mechanics of virus transfer between an infected and susceptible individual [23]. Using this model, we examined the transmission process’s sensitivity to the parameters associated with within-host or outside-host processes. The sensitivity analysis allowed for us to quantify the relative contribution of within-host and outside-host processes to the overall viral transmission.

2. Materials and Methods

2.1. Measuring Transmissibility

There are several ways to quantify the transmissibility of a new virus variant. This can be quantified as the effective reproduction number of the emerging new variant compared to that of the wild-type [24]. This comparison is often expressed in relative terms [20,25]. Another possibility is to describe transmissibility by estimating the secondary attack rate (SAR), defined as the probability of an infected person infecting a susceptible person [26,27]. Often, the SAR is stratified by the different venues in which people may be exposed to an infected person, e.g., public transportation, healthcare, and households [28,29]. Household settings were found to be associated with a high risk of infection by the SARS-CoV-2 virus [28,30]. This can be illustrated by the fact that, during the first wave of the COVID-19 pandemic in Israel, $65.8\%$ of the cases were infected at home [31]. The higher risk of infection can be attributed to the nature and duration of interactions between household members [32,33]. Household SAR (hSAR) is closely related to the reproduction number [34]. Several studies have characterized the increased transmissibility of a variant as the increase in hSAR relative to the wild-type. For example, the hSAR of the Alpha variant (lineage B.1.1.7, first detected in England [35]) was reported to increase by $55\%$ in Oslo, Norway, and the hSAR of the Epsilon variant (lineages B.1.427 and B.1.429, “California” or “West Coast” variants [35,36]) have increased by $40\%$ in San Francisco, United States [36,37]. Regarding the Delta variant (lineage B.2.627.2, first detected in India, [35]), the hSAR among Delta contacts in Singapore was $100\%$ higher than the hSAR among contacts exposed to other variants [38]. According to a report from the United Kingdom, in only a week since its first detection in late November 2021, the Omicron variant (lineage B.1.1.529, first detected in South Africa [39]) was sequenced in seven countries. The hSAR associated with this variant was reported to be $102\%$ higher than that of the Delta variant [40].

2.2. Mathematical Model

We utilized the mathematical model of Arav, et. al. [23] that describes the transmission of the SARS-CoV-2 virus between primary (infected) and secondary (susceptible) individuals in a household. The model accounts for the four transmission modes of the SARS-CoV-2 virus [41] (Figure 1): direct physical contact between people (Figure 1, mode 1); indirect physical contact via small, frequently touched objects (fomites), such as a doorknob or a faucet, (Figure 1, mode 2); touching environmental surfaces (large, less frequently touched surfaces such as table tops and other types of furniture), contaminated by droplets (Figure 1, mode 3); inhalation of droplet nuclei suspended in the air (Figure 1, mode 4).

Each mode of transmission comprises fast–discrete random events and slow–continuous events. The fast–discrete random events correspond to the actions of the individuals, such as touching their faces, touching each other, coughing, and talking. The slow–continuous events correspond to the degradation of virus copies on surfaces, hands, and in the air. The model equations calculate the kinetics of these two distinct dynamic regimes using a hybrid continuous, stochastic-jump framework [42]. In this framework, the actions of individuals (such as touching their faces or touching each other) are described as stochastic jump Poisson processes and represent behavior patterns that typically occur in the living room, kitchen, bath, and bedrooms. The degradation process of the virus on the hands, on surfaces, and in the air is described using continuum dynamics.

In the following, we describe the transfer of virus copies due to the individuals’ discrete actions Section 2.2.1), the kinetics virus concentration on surfaces, hands, and in the air Section 2.2.2), and finally, the probability of the secondary infection Section 2.2.3).

2.2.1. Virus Transmission Due to the Individuals’ Action

The transmission between the two individuals begins with the primary shedding of virus copies. Virus shedding occurs in the form of droplets and droplets-nuclei emitted when the primary sneezes, coughs, or talks or by depositing virus copies on the primary’s hands as the primary touches his face. The amount of viral shedding is proportional to the viral loads of the primary at that time ($L^P\left(t\right)$). During the incubation period, the viral load increases exponentially and reaches its maximal level ($L_{max}$) at its end [43]. The exchange of virus copies between the facial tissues to the fingers is given by:

$$\begin{array}{ccc}\hfill E_{mouthToFinger}^P\left(t\right)=& & {ϵ}_{skin-skin}\times V_{inoculation}\times L^P\left(t\right)\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill E_{FingerToMouth}^P\left(t\right)=& & {ϵ}_{skin-skin}\times C_{finger}^P\left(t\right)\times S_{finger}\hfill \end{array}$$

(2)

where $E_{MouthToFinger}^P\left(t\right)$ is the number of virus copies that pass from the facial membranes to the fingers of the primary at time t, $E_{FingerToMouth}^P\left(t\right)$ is the number of viruses that pass from the fingers to the facial membranes of the primary at time t, ${ϵ}_{skin-skin}$ is the fraction of viruses that transfer in skin-to-skin contact, $C_{finger}^P\left(t\right)$ is the concentration of the virus copies on the fingers of the primary in time t, $S_{finger}$ is the effective surface area of the finger, and $V_{inoculation}$ is the volume that is carried from the facial membranes by touch.

The virus is then transmitted from the hands of the primary to the hands of the secondary individual either directly (mode 1), indirectly via fomites (mode 2), or on environmental surfaces (mode 3).

Direct transmission occurs through physical contact with contaminated body parts. The number of virus copies transferred from the primary to the secondary in time t is given by,

$$E_{PhysicalContact}\left(t\right)={ϵ}_{skin-skin}\times \left(C_{finger}^P\left(t\right)-C_{finger}^S\left(t\right)\right)\times S_{finger}$$

(3)

where $C_{finger}^S\left(t\right)$ is the concentration on the fingers of the secondary in time t.

Indirect transmission occurs via fomites or environmental surfaces. The overall number of virus copies transferred from fomites to the hands following a single touch is given by:

$$E_{fomite}^I\left(t\right)=S_{finger}\times \left({ϵ}_{fomite-hand}\times C_{fomite}\left(t\right)-{ϵ}_{hand-fomite}\times C_{finger}^I\left(t\right)\right)$$

(4)

where $C_{fomite}$ is virus concentration on the fomite, ${ϵ}_{hand-fomite}$ and ${ϵ}_{fomite-hand}$ are the transfer coefficients from hand to fomite and vice versa, and I is the specific individual who touched the fomite, which can either be the primary or the secondary. The large droplets are emitted by the primary when he coughs, sneezes, or speaks contaminate the environmental surfaces. Each time the primary coughs, he creates a stain on the environmental surfaces. When touched, this stain would contaminate the hands of the individuals that touch it. The probability of touching a stain is proportional to its size ($S_{contaminatedArea}$) relative to the total surface area of the room. That is, the sum of the room’s surface area ($S_{room}$) and the furniture surface area ($S_{furniture}$). The number of virus copies that pass from the stain to the hands of the individual that touches it is given by

$$E_{surface}^I\left(t\right)=\sum _j^{N_{t,contaminatedArea}}{P}_j\times S_{finger}\times {ϵ}_{fomite-hand}\times C_{surface}\left(t\right)$$

(5)

where $P_j$ is the probability to touch a contaminated area, $N_{t,contaminatedArea}$ is the number of contaminated areas at time t, and $C_{surface}\left(t\right)$ is the virus concentration on environmental surfaces.

The inoculation of the secondary individual takes place after he touches his face (modes 1, 2, 3), as well as inhales virus-containing droplet nuclei (mode 4). The transfer rate from the secondary’s hands to his facial membranes is given by

$$E_{FingerToMouth}^S\left(t\right)={ϵ}_{skin-skin}\times C_{finger}^S\left(t\right)\times S_{finger}$$

(6)

where $E_{FingerToMouth}^s$ is the number of viruses that pass from the fingers to the facial membranes of the secondary, and $C_{finger}^S$, the concentration of the virus copies on the fingers of the secondary.

Following Nicas and Best [44], we assume that only a fraction ${ϵ}_{mouth}$ of the total number of virus copies deposited on the mouth reach the oro- and nasopharyngeal target sites. Another route of inoculation is the inhalation of the virus-containing droplet nuclei that were emitted (mode 4). Since only a fraction ${ϵ}_{breath}$ of the inhaled virus copies are deposited in the respiratory tracts [44], the contribution of this mode of transmission to the overall inoculation is given by:

$$E_{DropletNuclei}\left(t\right)={ϵ}_{breath}B·C_{air}\left(t\right)$$

(7)

where $E_{DropletNuclei}$ is exposure rate from the droplet nuclei, B is the breathing rate of the secondary individual, and $C_{air}\left(t\right)$ is the concentration of the droplet nuclei in the room in time t.

2.2.2. Kinetics of Virus Concentration

The probability that the secondary individual will become infected is calculated by his exposure to the virus copies he breathes and self-inoculates when he touches his face. To calculate the kinetics of exposure, we compute the kinetics of the concentration of virus copies on the hands of the primary, $C_{finger}^P\left(t\right)$, and the secondary, $C_{finger}^S\left(t\right)$, individuals, on fomites, $C_{fomite}\left(t\right)$, environmental surfaces, $C_{surface}\left(t\right)$, and in the air, $C_{air}\left(t\right)$.

The dynamics of the concentration on the fingers of the individuals is determined by inoculation, physical contact, touching fomites, environmental surfaces, hand washing, and by the degradation of the virus viability on the hands. Accordingly:

$$\begin{array}{ccc}\hfill \frac{d{C}_{finger}^I}{dt}=& & \sum _i^{{\mathcal{N}}_{inoculation}}\frac{E_{inoculation}^I\left(t_{i,inoculation}\right)}{S_{finger}}+\sum _i^{{\mathcal{N}}_{fomite}}\frac{E_{fomite}^I\left(t_{i,fomite}\right)}{S_{finger}}\hfill \\ & & +\sum _i^{{\mathcal{N}}_{surface}}\frac{E_{surface}^I\left(t_{i,surface}\right)}{S_{finger}}\hfill \\ & & +\delta \sum _i^{{\mathcal{N}}_{PhysicalContact}}\frac{E_{PhysicalContact}\left(t_{i,PhysicalContact}\right)}{S_{finger}}\hfill \\ & & -\sum _i^{{\mathcal{N}}_{washing}}{E}_{washing}^I\left(t_{i,washing}\right)-{\alpha }_{hand}\times C_{finger}^I\hfill \end{array}$$

(8)

where I is either the primary or the secondary individual. The first 5 terms on the right hand side of Equation (8) describe the transfer of virus copies as a result of discrete events and the last term corresponds to the degradation of the virus viability on the hands where ${\alpha }_{hand}$ is the degradation rate constant. The $\delta$ equals 1 for the primary individual and $-1$ for the secondary individual. The discrete event times ($t_{i,X}$ with $i=1‥{\mathcal{N}}_{\mathcal{X}}$ where X is inoculation, fomite, surface, washing and physical contact) are determined from a Poisson distribution with rate constant ${\tau }_X$.

The concentration of virus copies on fomites, $C_{fomite}$, is determined by touching fomites and the degradation of virus viability:

$$\begin{array}{ccc}\hfill \frac{d{C}_{fomite}}{dt}=& & -\sum _i^{{\mathcal{N}}_{fomite}^P}\frac{E_{i,fomite}^P\left(t_{i,fomite}\right)}{S_{fomite}}-\sum _i^{{\mathcal{N}}_{fomite}^S}\frac{E_{i,fomite}^S\left(t_{i,fomite}\right)}{S_{fomite}}\hfill \\ & & -{\alpha }_{fomite}\times C_{fomite}\hfill \end{array}$$

(9)

where ${\alpha }_{fomite}$ is the degradation rate on the fomite.

The concentration of the virus copies in the air is calculated from the emission of droplet-nuclei during the coughing, sneezing, and talking activities of the primary individual and the ventilation process:

$$\begin{array}{ccc}\hfill \frac{d{C}_{air}}{dt}=& & \sum _i^{{\mathcal{N}}_{cough}^P}\frac{V_{cough,aerosol}\times L^P\left(t_{i,cough}\right)}{V_{room}}+\sum _i^{{\mathcal{N}}_{sneeze}^P}\frac{V_{sneeze,aerosol}\times L^P\left(t_{i,sneeze}\right)}{V_{room}}\hfill \\ & & \sum _i^{{\mathcal{N}}_{talk}^P}\frac{V_{talk,aerosol}\times L^P\left(t_{i,talk}\right)}{V_{room}}-\beta \times C_{air}\hfill \end{array}$$

(10)

where $C_{air}$ is the virus concentration in the air and $L^P$ is the viral load of the primary individual, ${\mathcal{N}}_{cough}^P$, ${\mathcal{N}}_{sneeze}^P$ and ${\mathcal{N}}_{talk}^P$ are the number of times that the primary coughs, sneezes, and talks (respectively). These occur at time $t_{i,cough},t_{i,sneeze}$ and $t_{i,talk}$, and are determined by a Poisson distribution with time constants of ${\tau }_{cough},{\tau }_{sneeze}$, and ${\tau }_{talk}$, respectively. The dynamics of the concentration in the k-th-contaminated area on an environmental surface that was created by a cough, sneeze or during a talk activity at time $T_k$ is determined by the following equation:

$$\begin{array}{ccc}\hfill \frac{d{C}_{surface}^k}{dt}=& & -\sum _i^{{\mathcal{N}}_{surface}^P}P\times \frac{E_{i,surface}^{k,P}\left(t_{i,surface}\right)}{S_{contaminatedArea}}-\sum _i^{{\mathcal{N}}_{surface}^S}P\times \frac{E_{i,surface}^{k,S}\left(t_{i,surface}\right)}{S_{contaminatedArea}}\hfill \\ & & -{\alpha }_{surface}\times C_{surface}^k\hfill \end{array}$$

(11)

where $C_{surface}^k$ is the concentration on the k-th surface, $P=\frac{S_{contaminatedArea}}{S_{room}+S_{furniture}}$ is the probability to touch a contaminated area and ${\alpha }_{surface}$ is the degradation rate in this surface.

The k-th-contaminated area was created at time $T_k$ by a cough, sneeze or during speech. The initial concentration in the contaminated area is the product of the viral load of the primary individual at time $T_k$ with the volume of droplets that are larger than 100 µm ($V_{X,droplets}$, where X is cough, sneeze or talk).

The exposure of the secondary individual is calculated within the exposure time interval and calculated as the sum of exposure in the contact and droplet nuclei mode of transmission. Following Nicas and Best [44], we assume that only a fraction ${ϵ}_{mouth}$ of the total number of virus copies deposited on the mouth reach oro- and nasopharyngeal target sites, and only a fraction ${ϵ}_{breath}$ of the inhaled viruses deposited in the respiratory tracts. Therefore, the exposure is given by the equation:

$$\begin{array}{ccc}\hfill \frac{d{D}_{secondary}^i}{dt}=& & {ϵ}_{mouth}\sum _i{E}_{FingerToMouth}^{secondary}\left(t_{i,inoculation}\right)+E_{DropletNuclei}\left(t\right)\hfill \end{array}$$

(12)

where $D^i$ is the exposure obtained when solving the equation between the two exposure events that take place at times $T_{exposure}^i$ and $T_{exposure}^{i+1}$ and B is the breathing rate.

2.2.3. Probability of Infection

The probability that a secondary individual will be infected depends on the number of virus copies to which he was exposed and the probability of getting infected by a single exposure dose (measured by plaque forming unit, PFU, or viral copies) within a particular time interval. We divide the simulation into periods and calculate the overall exposure for each period. The probability that a susceptible individual will become infected, $P_{infected}$, was inferred from the dose–response curve first reported for SARS-CoV-1 and was later found to agree well with exposure data to SARS-CoV-2 [45,46].

$$P_{infected}=1-e^{-\frac{D_j}{\lambda k}}$$

(13)

where $D_j$ is the total exposure of the individual at time period j (see Equation (12) below), k is the reciprocal in of a single plaque-forming units (PFU) to get infected, $\lambda$ is a conversion factor that expresses the number of virus copies in a single PFU. The length of each time period was determined from a Poisson process, using a timescale of 7 h that represents the life cycle of the virus in the body, inferred from the life cycle of the SARS-COV-1 [47].

2.2.4. Model’s Results and Parameters

A single realization of the model, r, begins when the primary becomes infected and begins an incubation period whose duration, $I_r$, is drawn from a log-normal distribution with a mean of 5 days and standard deviation (SD) of $0.45$ days [48]. The viral load of the primary increases exponentially with time [49], reaching a maximal level at the end of the incubation period [43]. The simulation ends when the secondary individual is infected or after $I_r$ days. Two variables. The results of the realization r, ${\mathsf{\Lambda }}_r$ and $t_{\mathsf{\Lambda }}$.

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_r=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{t}_{{\mathsf{\Lambda }}_r}\ge 0\hfill \\ 0\hfill & otherwise\hfill \end{array}\right.\end{array}$$

(14)

where $t_{{\mathsf{\Lambda }}_r}$ is the time in which the secondary became infected and ${\mathsf{\Lambda }}_r$ is a binary variable, which signifies whether the secondary individual became infected during the realization r.

Since the actions of the individuals are represented as a stochastic process, a Monte Carlo simulation was performed, in which multiple realizations were computed to obtain the appropriate ensemble statistics, hSAR, and SI, and the serial interval (SI) distribution. The household secondary attack, hSAR, is defined as follows:

$$\begin{array}{ccc}\hfill hSAR=& & \frac{{\sum }_{i=1}^R{\mathsf{\Lambda }}_r}{R}\hfill \end{array}$$

(15)

where R is the total number of realizations.

The SI is the distribution of the time period between the symptoms’ onset in the primary and the secondary. This distribution is associated with the reproductive number and other fundamental variables in epidemiological models [50,51]. According to the many studies, this distribution is considered to follow the gamma distribution, characterized by its mean and standard deviation. Estimates for these values for the wild-type variant range between 4.03 and 6.3 days and a standard deviation of 3–4.2 days [30,52,53,54]. In this study, the SI is the empirical distribution of $t_{\mathsf{\Lambda }}$.

The model parameters were taken from the available literature and are described in the Appendix A.

2.3. Modelling Effect of Within-Host and Outside-Host Processes on Transmissibility

Each transmission mode consists of within-host and outside-host processes. The outside-host processes are the degradation of the virus copies in the air and on hands, fomites, and environmental surfaces. The results of the within-host dynamics are the contagiousness of the primary and infectivity of the secondary. The model describes the contagiousness using the primary’s maximal viral load since this parameter is proportional to the extent of his viral shedding [55,56,57]. The infectivity is represented in the model by the probability of the secondary individual becoming infected by a single exposure dose and corresponds to $\frac{1}{k}$[46]. This is because an increase in the likelihood of infection is a decreased infectious dose, k, meaning that a smaller number of virus copies is required to initiate an infection [4]. It is noteworthy that this increased likelihood of infection can also be considered as the increased susceptibility of the secondary [58].

The mathematical model outlined above (Figure 1) was used to examine the sensitivity of within-host and outside-host processes on transmissibility. We used the model to calculate the changes in hSAR and the SI of different hypothetical variants with enhanced outside-host or within-host behavior.

Two kinds of variants with enhanced outside-host behavior were examined. Variants of the first kind had a higher stability (or slower degradation) on surfaces, and the second variants had a higher stability within the droplet nuclei suspended in the air. Specifically, the effect of higher stability on surfaces was examined by decreasing the degradation rate on furniture (${\alpha }_{furniture}$), fomites (${\alpha }_{fomite}$), and on the hands of the primary and the secondary individuals (${\alpha }_{hand}$). The effect of higher stability in the air was examined by decreasing the degradation rate of the virus in the aerosol (${\alpha }_{air}$). The benchmark values that represent the wild-type were ${\alpha }_{air}$, ${\alpha }_{furniture}$, and ${\alpha }_{fomite}$ of 1 h${}^{-1}$, $6$ h${}^{-1}$, and 6 h${}^{-1}$, respectively [59]. To quantify the effect of these two kinds of variants relative to the wild-type, the benchmark degradation rates were divided by values ranging from 2 up to 8. This range corresponds with reports on a 10-fold range of aerosols containing different strains of foot-and-mouth disease virus and the broad range of stability rates of different Influenza A strains on some surfaces [60,61]. It should also be noted that ${\alpha }_{furniture}$ and ${\alpha }_{fomite}$ were set to equal values, based on the fact that the difference between environmental and fomite surfaces lies in their area and the rate that each type of surface is touched, but not in the degradation rate of the virus copies that are deposited on them [23].

Regarding the within-host behavior, two kinds of variants were examined. The first kind was a more contagious variant, characterized by a higher maximal viral load ($L_{max}$). The second kind was a more infectious variant, characterized by a higher probability of infecting with a single exposure dose (that is, a higher reciprocal of the dose–response coefficient, k). The benchmark value of $L_{max}$ that represents the wild-type SARS-CoV-2 virus was 2·${10}^8$ copies [62] and a k of 410 PFU, which translates to $1.2$·${10}^5$ viral copies [45,46]. These two types were examined by increasing the wild-type benchmark values from $10\%$ to $100\%$. That range is well within the reported range of the minimal dose required to cause an infection or antibody response by different strains of foot-and-mouth disease virus and adenovirus that spans over a two-fold range to an order of magnitude, respectively [63,64].

2.4. Numerical Method

To solve the model’s Equations (8) to (12) we chose to use the jump-adapted approximation proposed by Casella et al. [65]. For each realization, we determined the discrete event times ($t_{i,X}$ with $i=1\cdots {\mathcal{N}}_X$ where X is inoculation, fomite, surface, washing, coughing, physical contact, sneezing and talking) are determined from a Poisson distribution with rate constant ${\tau }_X$. Then, the events were combined and sorted in ascending order to obtain a set of discrete times. That is, $\mathcal{T}=\{T_1,T_2,T_3,\cdots ,T_{{\mathcal{N}}_{total}}\}$, where each $T_i$ is assigned to a one $t_{i,X}$ event and ${\mathcal{N}}_{total}$ is the total number of discrete events.

The simulation is then implicitly solved between the times $T_i$ and $T_{i+1}$ with a time step of 6 seconds. At each $T_i$, we calculate the number of virus copies transferred for the corresponding event that takes place at this point.

The model was implemented in python version 3.6.5 and is available online (https://www.python.org/downloads/, accessed on 18 July 2022), as detailed in the data availability statement.

3. Results

3.1. Convergence Analysis

The model’s Monte Carlo simulation consists of multiple realizations, from which the hSAR and SI are calculated. To determine the value of R, the total number of realizations required for the convergence of the results, a preliminary step was performed to examine the convergence of the model’s simulation. A simulation with R of 10,000 was conducted. Then, we performed a Monte Carlo simulation with different values of R. For each simulation, the relative error was calculated:

$$relativeerror\left(R\right)=\frac{\left|hSAR\right(R)-hSAR(\mathrm{10,000}\left)\right|}{hSAR\left(\mathrm{10,000}\right)}$$

(16)

where $hSAR\left(R\right)$ is the hSAR calculated by Equation (15) for R realizations. The relative error of hSAR began with values of more than $40\%$ for R values of a few hundred realizations. It then decreases with the increase in R, where the rate of decrease results in a convex graph. At R of 4000 the graph becomes quasi-stationary, where larger R values all result in relative errors of around $1\%$ (Figure 2).

3.2. Reaction to Changes in Within-Host or Outside-Host Processes

3.2.1. Effect on Household Secondary Attack Rate

We begin by examining the effect of outside-host processes on the hSAR, as was delineated in Section 2.3. Figure 3A shows the increase in the hSAR of variants with a higher stability on surfaces. A relative decrease of up to a factor of 2 in the degradation rates of either kind of variant leads to a relative increase in the hSAR, but only up to $10\%$. However, a further decrease in the degradation rates does not lead to a further rise in the hSAR. Figure 3B shows the increase in teh hSAR of variants with higher stability in the droplet nuclei suspended in the air. A higher virus stability in droplet nuclei yields a similar rise in the hSAR, resulting from higher surface stability. Specifically, the hSAR increases by, at most, $10\%$. This relative increase is far from the reported relative hSAR increase of $40\%$ associated with the Epsilon variant, let alone from the even more transmissible variants, Alpha, Delta, and Omicron [36,37,38,40].

Next, we examine the effect of within-host processes on the hSAR (see Section 2.3 for details). Figure 3C shows the effect of genetic variations that increase the contagiousness on the hSAR. As seen, the relative increase in contagiousness is linearly associated with the relative increase in the hSAR. Specifically, the fitted linear model that describes this effect yields a slope of $0.68$ and an intercept of $0.342$ ($R^2=0.997$). Hence, an increase of $56\%$, $78\%$, and $144\%$ in the maximal viral load, relative to the wild-type, could account for the reported increase in the transmissibility of the Epsilon, Alpha, and Delta variants ($40\%$, $55\%$, and $100\%$, respectively) [36,37,38]. Regarding the Omicron variant, to account for its reported $102\%$ increase in the hSAR, the maximal viral load should increase by $129\%$ compared to the Delta variant [38,40]. Figure 3D shows the effect variants with higher infectivity on the hSAR. Similarly, a linear relationship was also found between the infectivity and the hSAR, with a slope of $0.79$ and an intercept of $0.211$ ($R^2=0.998$). Hence, the probability of infection by a single exposure dose should increase by $51\%$, $69\%$, $126\%$, and $129\%$, respectively, to account for the higher transmissibility reported for the Epsilon, Alpha, Delta, and Omicron variants.

3.2.2. Effect on the Serial Interval

The SI of variants with enhanced outside-host and within-processes resulted in an SI that is similar to the wild-type. For example, a relative increase of $700\%$ in the stability of droplet nuclei suspended in the air (Figure 4C) or a $600\%$ increase in the stability of virus copies on surfaces (Figure 4D) resulted in SI that is close to the SI of the wild-type. Both are within the boundaries of the observed SI (Figure 4, shaded area). Furthermore, a relative increase of $150\%$ in the contagiousness (Figure 4A) and $100\%$ in the infectivity led to a similar outcome.

4. Discussion

In this work, we examined whether the enhanced transmissibility of the recently emerging SARS-CoV-2 variants could result from within-host or outside-host processes. This was achieved by computing the increase in the hSAR of several hypothetical variants with Monte-Carlo simulations of a well-established mathematical model for the transmission of SARS-CoV-2 in households.

In contrast to the limited effect of outside-host processes, we have found that variants with enhanced within-host processes linearly increase hSAR (coefficient of 0.7–0.8, Figure 3C,D). Specifically, it is expected that a variant with higher contagiousness or higher infectivity (or the combination of the two) will increase its transmissibility to the observed range of the Epsilon, Alpha, Delta, and Omicron variants.

It is interesting to associate our finds with a couple of recent studies [9,36]. These have reported that the viral load associated with the Alpha and Epsilon variants was not significantly different from the viral load of other circulating strains. Therefore, these variants are not characterized by higher contagiousness. Hence, we conclude that higher infectivity accounts for the increased transmissibility of these SARS-CoV-2 variants. While the exact mechanism by which genetic variations increase the infectivity is not currently known, Nelson, et al. [12] and Gan, et al. [66] have shown that certain variations in the spike RBD increase the affinity to the human ACE2 and thereby promote cell entry. It seems reasonable that the manifestation of this affinity is enhanced infectivity. Contrary to the viral load results regarding the Alpha and Epsilon variants, the viral load associated with the Delta variant was at least ten times higher than the historical variants [8,67]. Therefore, in the case of the Delta variant, it seems that it is driven by higher contagiousness. Unlike these variants, which by the time scales of the COVID-19 pandemic have circulated for a relatively long time, the spread of the Omicron variant is an ongoing event. Therefore, the available literature is limited, and most has been published on pre-print servers and should be subjected to peer review. One of these studies analyzed the time series of the number of people that were reinfected by SARS-CoV-2 in South Africa. It was found that the emergence of the Omicron variant in South Africa temporally coincides with a substantial and ongoing increase in the daily number of new reinfections. This suggests that this variant is associated with an ability to evade the immunity induced by a prior infection [39]. Another study used a variety of computational tools to examine the Delta and Omicron variants. It found that, compared to the binding affinity for human ACE2, which is facilitated by the Delta variant’s RBD, the numerous mutations in the Omicron variant’s RBD promote a higher binding affinity [11]. As both immune evasion and higher ACE2 affinity can be described as an increased probability of getting infected by a single virus copy, the current information leads to the conclusion that the higher transmissibility of the Omicron variant can be attributed to its enhanced infectivity.

The modeled SI results were found to be rather similar to that of the wild-type. These results apply to outside-host processes as well as within-processes. Some authors mention that the specific reason for the fast rise in the number of cases infected by the Omicron variant compared to the Delta variant was unclear. Possible reasons were an increase in transmissibility (which can be expressed by hSAR, as described in Section 2.1), shortening of the SI, or a combination of these [68,69]. It can be argued, with some degree of caution, that the SI results of this study provide evidence that can settle the debate towards the increase in transmissibility explanation. However, it should be noted that several studies have indicated that the mean SI of the Omicron variant is shorter compared to the Delta variant. This shortening ranges between 0.3 and 0.6 days [68,69,70,71].

The results of this study emphasize the importance of within-host processes in the spread of highly transmissible SARS-CoV-2 variants. Therefore, a prospect for future research could be expanding the model used in this study to provide a further, detailed description of within-host processes. Several models follow the spread of SARS-CoV-2 infection within a host. These models follow the dynamics of one or several variables or species. The species that all these models describe is the number of virus copies within the host. Models that also describe other species may recount the numbers of the host’s target and immune cells in various stages, RNA genome, and cytokines stimulating inflammatory responses [72,73]. Given the introduction of highly transmissible SARS-CoV-2 variants, a significant finding was the phenomenon of within-host diversity. This means that, within a single infected individual, numerous viral mutations may develop. However, most of these mutations are not transferred to other hosts and, therefore, are lost due to the genetic bottleneck [74]. This effect is defined as the number of virus copies that genetically contribute to the infection process. To incorporate this effect in a model, it is necessary to estimate the probability of infection caused by the number of virus particles contained within the size of the genetic bottleneck [75]. Such an estimate may be provided by virus bioinformatics tools, some of which were recently developed to address the needs of SARS-CoV-2 research. The Haploflow graph-based assembler, which allows for detection and full-strain reconstruction, is notable. Thus, it is appropriate to study SARS-CoV-2 infection scenarios that involve multiple viral strains [76].