The transmission of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the agent that causes coronavirus disease 2019 (COVID- 19), is dependent on the amount of infectious particles present in the environment surrounding the susceptible host, and/or on the proximity between the susceptible and infectious host. Experimental studies that use real-time reverse transcription-polymerase chain reaction (PCR) assays have reported the presence of SARS-CoV-2 RNA in contaminated environmental surfaces [1][2][3][4] and in aerosols [5][6][7][8]. Moreover, in vivo cell culture assays, have shown that particles released into the environment are replication-competent [9] and they can stay infectious in the aerosols for up to three hours [10].
Once an infection is established, the dynamics of SARS-CoV-2 in the upper respiratory tract (URT) is dictated by the interplay between the virus fitness and host immune responses. The total RNA to infectious virus ratio changes over the course of the infection, varying between 10 3 : 1 and 10 6 : 1 RNA to plaque forming units (PFU) [11,12]. Determining the within-host mechanistic interactions responsible for the temporal changes in infectious to non-infectious viral dynamics is important for guiding interventions.
Over the last two years, within-host mathematical models developed for influenza and other respiratory infections have been modified for SARS-CoV-2 infections [13][14][15][16][17][18]. These models divided total viral titers into infectious and non-infectious particles [19,20], fitted their sum to total RNA values measured by PCR (used as a proxy for total virus load) and used the results to determine the mechanisms of viral expansion and loss. The models were subsequently used to provide insights into the relative fitness of variants, types of drug interventions [13][14][15], the relationship between individual infection and population transmission, and the effect of this relationship on testing and vaccine strategies [21,22]. In most studies, the total RNA to infectious virus ratio in URT is assumed constant over time. However, Goyal et al. [23] and Ke et al. [13] proposed a nonlinear correspondence between infectious virus and total RNA, and found that a density dependent function best describes their relationship [13].
While there is a reasonable understanding of the mechanistic dynamics modulating SARS CoV-2 infection in the URT, the temporal shedding into the environment of viral RNA and infectious virus has not been explored. Here, we expand a within-host mathematical model to include the dynamics of viral RNA and infectious virus titers in both URT and aerosols. We validate the models against two URT and one aerosol inoculation study in golden hamsters [24,25] and use the models to determine the relationship between infectious virus and total RNA in both environments. Lastly, we investigate the link between the infectious virus shed in the environment and the probability of a nearby host getting infected. The results can guide interventions.
We use previously published temporal SARS-CoV-2 RNA and infectious virus titer data from two inoculation studies in golden Syrian hamsters:
• Sia et al. study [25]: six donor hamsters (all male) were inoculated intranasally with 8 × 10 4 tissue culture ineffective dose (TCID 50 ) of SARS-CoV-2. At 24 hours after inoculation, each donor was transferred to a new cage and co-housed with one naive male hamster. Viral RNA and infectious virus titers were collected every other day (in RNA/ml and TCID 50 /ml) for the first 14 days in both donors and contacts and their weight changes were monitored daily. We will refer to these two groups as donors and contacts.
• Hawks et al. study [24]: eight hamsters (4 males and 4 females) were inoculated intranasally with 10 5 PFU (1.4 × 10 5 TCID 50 ) of SARS-CoV-2. Viral RNA and infectious virus titers were collected daily from nasal washes (in RNA/wash and PFU/wash) and from the exhaled breath (aerosols) (in RNA/hour and PFU/hour) for the first five days and then again at day 10. Since the nasal wash and air samples were collected in 100 μl and 400 μl, we rescaled the RNA/wash data to RNA/ml by multiplying it with 10 and 2.5, respectively. Moreover, since 1 TCID 50 = 0.7 PFU, we rescaled the PFU/wash infectious virus titers in the nasal washes and in the air to TCID 50 /ml by multiplying them with 10/0.7 and 2.5/0.7, respectively. Hamsters were weighed daily and the study was terminated when clinical signs of illness were observed. We will refer to these two groups as males and females.
We model the interaction between target epithelial cells
with initial conditions T(0) = T 0 , E(0) = 0, I(0) = 0, V u (0) = V 0 , V a (0) = 0. Moreover, we model total viral RNA in the upper respiratory tract R u and in the air R a as follows. Since the measured viral RNA may not be a good representation of actual virions (infectious or noninfectious) but also of naked viral RNA from virus neutralization or infected cells death, we assume viral RNA is produced at rate ρ 1 per infected cell, and rate ρ 2 per infectious virus. Rate ρ 1 incorporates both newly synthesized genomic RNA that can be used for replication or transcription, p, and naked RNA released from dead infected cells, q, so that ρ 1 = p + q. Rate ρ 2 represents neutralized infectious virus. Genomic RNA particles are removed at rate d + c (where d is the degradation rate and c is the immune clearance rate, as before), while naked RNA is lost at rate d (with the immune system not being involved in its removal). Lastly, URT total viral RNA is released into the environment at rate ϕ 2 where it is lost at rate d. The equations describing total viral RNA dynamics over time are given by
We assume an initial target cell population in the upper respiratory tract T(0) = 10 7 epithelial cells/ml, as in influenza [26], and no exposed and infected cells E(0) = I(0) = 0 epithelial cells/ ml. The initial infectious virus is given by the inoculum titer, V u (0) = 8 × 10 4 TCID 50 for donors and V u (0) = 1.4 × 10 5 TCID 50 for males and females. For contacts, who are infected rather than inoculated, V u (0) is unknown. We investigate two scenarios, (i) infection starts with a single infectious virus, V u (0) = 1 TCID 50 , and (ii) infection starts with more than one infectious virus V u (0)>1 TCID 50 , with the exact number being estimated through data fitting. No viral RNA is present in URT at the time of inoculation, V a (0) = 0 TCID 50 , and neither infectious virus nor viral RNA are present in the air at the time of inoculation, R u (0) = R a (0) = 0 TCID 50 . We assume that the infectivity rate is β = 5 × 10 -6 mL/(TCID 50 x day), infectious virus removal rate is c + d = 10 per day [13], the eclipse rate is k = 4 per day [13], the RNA degradation rate is d = 1 per day, the RNA release due to neutralization is ρ 2 = c = 9 per day, and the enhanced inactivation due to the elements is d 1 = 1. Lastly, since RNA and infectious virus from exhaled breath are only collected in Hawks et al. [24], we assume ϕ 1 = ϕ 2 = 0 in models (1) and (2) when applied to donors and contacts and {ϕ 1 , ϕ 2 } 6 ¼ 0 when applied to males and females.
The weights of males and females vary from 56 to 76.1 grams (with an average of 70 grams) at inoculation. After challenge, each subject was weighed, as a way to monitor them for clinical signs of illness, and daily percent weight changes were reported as changes from the baseline weight w 0 = 1, 0 � w i � 1 for days i 2 S = {1, 2, 3, 4, 5, 10} post infection. To account for all intermediate time points, we fitted a 4-degree polynomial to the weight data (smallest degree polynomial that gave a residual sum of squares <10 -3 ), and the resulting percent weight functions are shown in Fig 2A.
The initial weights of donors and contacts are unknown. After inoculation, each subject was weighed and percent weight changes were reported as changes from the baseline weight w 0 = 1, 0 � w i � 1 for days i 2 S = {2, 4, 6, 8, 10, 12, 14} post infection. As before, we fitted a 4-degree polynomial to the weight data, and the resulting percent weight functions are shown in Fig 2B and2C for the donors and contacts, respectively.
To include the weight variability into the model, we assume that weight changes result in adjusted available epithelial cells [27]. At day t target, exposed and infected cells become and
Moreover, the initial conditions become
0 where ξ = {1.06, 1.09, 0.97, 0.81, 1.01, 0.92, 1.094, 1.04} is scaling from the average initial weight of 70 grams for males and females and ξ = 1 for donors and contacts [25].
The basic reproduction number (or basic reproductive ratio) is defined as the number of infected cells (or virus particles) that are produced by one infected cell (or virus particle) when the virus is introduced into a population of uninfected target cells T 1 (0) = ξT 0 . It is given by
with ϕ 1 = 0 for donors and contacts.
Using models ( 3) and ( 4) we estimate parameters param 1 w ¼ fd; p; � 1 ; r 1 ; � 2 g, by minimizing the functional
where
; and S = {1, . . ., 5, 10} days post infection.
Using model (3) we estimate parameters param 2 w ¼ fd; p; rg for donors and either param 2 w ¼ fd; p; rg with V u (0) = 1 TCID 50 or param 3 w ¼ fd; p; r 1 ; V u ð0Þg for contacts by minimizing the functional
; S = {2, 4, 6, 8, 10, 12, 14} days post infection for donors, S = {1, 3, 5, 7, 9, 11, 13} days post infection for contacts, and i = {2, 3}. We only consider the first data point at or below limit of detection. We use the 'nlinfit' algorithm in matlab and the resulting estimates for parameter means and 95% confidence intervals are given in Tables 1 and2. The theoretical solutions of V u and R u versus male and female data are shown in
We assume that the log10 URT infectious virus titer, ν = log 10 V u (measured in log10 TCID 50 per mL) can be modeled as a density dependent function of the total viral RNA, R u (measured Individual estimates (mean and 95% confidence intervals) from simultaneously fitting V u and R u given by model Eq (3) to URT infectious virus and RNA data in donors and contacts from Sia et al. [25].
https://doi.org/10.1371/journal.pcbi.1009997.t002
in RNA per mL), as in Ke et al. [13] and Goyal et al. [23])
To avoid overfitting, we fix V m = 8.5 and estimate the remaining parameters param 4 w ¼ fh; K m g by fitting Eq (6) to the population (ν, R u ) data from donors, contacts, males URT, and males aerosol. We excluded the female group due to limited infectious virus titers above the limit of detection. All data at and below limit of detection is treated as censored during data fitting. We maximize the functional,
using the 'nlinfit' algorithm in The resulting population parameter means 95% confidence intervals are given Table 3 and the population fits for log10 infectious virus, ν = log 10 V u , versus total viral RNA, R u , for each group are shown in Fig 5A -5D.
To study the kinetics of infectious virus titers and total viral RNA in the upper respiratory tract we used a dynamical within-host target cell limitation model developed for other acute respiratory infections (Eq 3), which was normalized to include the temporal change in the subject's weight (see Materials and methods for a full description). We estimated unknown biological parameters by fitting the model to infectious viral titer and total RNA data of four groups (males, females, donors, contacts) from two inoculation studies [24,25] (see Materials and methods). The resulting dynamics are in good agreement with the data kinetics in all four groups, with some inter group differences in the predicted outcomes. Infectious virus titers in the male and donors groups, who were challenged with high viral dose, expand to reach a peak 12-24 hours after inoculation and decline below limit of detection by 5-10 days post inoculation (p. In the hamsters of the contact group, who were infected (rather than inoculated) the inoculum dose is unknown. We assumed that infection is started by a single virion V u (0) = 1 TICD50 and found that the infectious virus titers have a delayed expansion compared to the other three groups, peaking 1.5-2.4 days post infection and decaying below limit of detection 7. We find consistent estimates between the groups for the infected cell death rate δ, with ranges between 1.5-5.1 per day, corresponding to infected cells life-spans of 4-16 hours. The production rate of infectious virus is highly variable among groups. In the female group we fixed the production rate to p = 0.1 TCID 50 /infected cell. This was significantly lower than the estimated production rates in the other groups, with contacts, donors and males average production rates of p = 20.4, p = 21.7 and p = 45 TCID 50 /infected cell, respectively. Moreover, the RNA production rate is variable among the two studies but similar among the groups within 1 and2. https://doi.org/10.1371/journal.pcbi.1009997.g006 each study. In particular, ρ 1 = 161 RNA/infected cells and ρ 1 = 494 RNA/infected cells in males and females, and two orders of magnitude lower than in donors and contacts, ρ 1 = 9 × 10 4 RNA/infected cells, respectively. The differences between the studies can be noted in the higher total RNA in donors and contacts (Fig 4 blue curves), 2.5 order of magnitude higher than in males and females (Fig 3 blue curves). Overall infectious virus titers (and consequently theoretical predictions) are two fold higher in donors and contacts compared to males and females, a difference we attribute to the experimental setup in [25] versus [24].
In order to determine the relationship between the amount of infectious virus titer (and total viral RNA) shedding into the environment over time and host-viral dynamics in the URT, we added two viral compartments to the within-host model and normalized them to account for subject's weight variability. The resulting within-host and aerosols model is given by systems Eqs 3 and 4 (see Materials and methods for full derivation). We estimated viral shedding parameters by the models against temporal aerosol data (infectious and total RNA) from the males and females (see Materials and methods). The resulting dynamics are in good agreement with the aerosol kinetics in both groups. In particular, the models predict that both infectious virus titers and total viral RNA get released into the air immediately after inoculation. For the males, the shedded infectious virus titers peak 21-37 hours after inoculation, and decay below limit of detection 3.7-5.6 days after inoculation (Fig 3A right panels, red lines). Total viral RNA peak 1.5-2.5 days post inoculation, and persist above limit of detection for the duration of the experiment (Fig 3A right panels, blue lines). For the females, the shedded infectious virus titers are maximal 3-5 hours after inoculation and decay below limit of detection 2.5-4.1 days later, faster than in the male group (Fig 3B right panels, red lines). Females total viral RNA kinetics, however, are similar to those in male, peaking 1.5-1.8 days after inoculation and persisting above limit of detection for the duration of the experiment (Fig 3B right panels, blue lines). We observe sex-specific differences in the shedding rates, with females infectious virus shedding rate ϕ 1 = 6.6 × 10 -3 being 132-times higher than that of the males, ϕ 1 = 0.5 × 10 -4 per day. The RNA shedding rates are similar among the sexes, ϕ 2 = 7.3 × 10 -5 and 11 × 10 -5 per day for males and females, respectively.
For each group, we estimated the within-host basic reproductive number R 0 , which averaged at 61 for the male group, 0.3 for the female group, 44 for the donor group, and 41 for the contact group.
Throughout the course of an infection hamsters shed both infectious and non-infectious virus into the air. Transmission to a close contact occurs when the infectious viruses reach the recipient and establish an infection. We wanted to determine whether aerosol data is a good predictor for the number of infectious virions that jump start such an infection in a close contact. While we know the exact inoculum value for the hamsters in the donor group, we do not know the inoculum value for the hamsters in the contact group. We first assumed that one infectious virus V u (0) = V 0 = 1 TCID 50 is sufficient to start the infection in the contact group. This lead to delayed viral expansion compared to the donor group, but similar clearance time and R 0 values. Next, we investigated how these results change if we assume that a larger number of infectious virions are needed to start the infection in the contact group. We included the infectious inoculum V u (0) = V 0 as an unknown parameter, and estimated it (together with parameters parm We found that the inoculum varies over a range between 40-1820 TCID 50 /ml among the six contacts. Since the contacts were co-housed with the infected donors at day 1 post inoculation, we compared the estimated V u (0) with the amount of infectious virus found in the aerosols of males and females. Our model predicted that the aerosol values for the infectious virus titers at day one post inoculation, V a (1), ranged between 96-222 TCID 50 /ml in males and between 7-52 TCID 50 /ml in females (see Fig B in S1 Text). When accounting for the fact that measured infectious virus in the Sia et al. study [25] is two-folds higher than the measured infectious virus in the Hawks et al. study [24], we find that the aerosol measurements are a good proxy for infectious inoculum. We have, however, found large confidence intervals for the estimated V u (0) in all contact hamsters (Table A in S1 Text), suggesting that more information is needed to determine the role of the infectious inoculum in the observed outcomes.
Detection of viral RNA by PCR testing from URT samples is the gold standard for COVID-19 diagnosis and is used to instate and discontinue control precautions, such as isolation and quarantine. However, there is no clear correlation between detection of viral RNA and detection of infectious virus titers. We use the data in the two hamster infection studies [24,25] and our model predictions to investigate: (i) the dependence of infectious viral levels on the RNA levels, and (ii) the dynamics over time of the viral RNA to infectious virus ratio. We exclude the female group from these analyses due to limited amount of infectious virus data above the limit of detection.
To quantify the dependence of infectious viral levels on the viral RNA levels we combined all population (R u , V u ) pairs. This resulted in 14 donor, 19 contact, 15 male URT, and 12 male aerosol pairs above the limit of detection and 24 donor, 17 donor 6 male URT, and 9 male aerosol pairs below the limit of detection. We assumed that the two quantities can be described by a density dependent function
(Eq 6 in Materials and Methods) as in [13,23]. When we fitted this function to the population (R u , V u ) data in the four groups (taking into consideration the presence of censored data, see Materials and methods for full explanation), we found that the level of infectious viruses increases sub-linearly with increases in viral RNA, with estimated exponent h ranging between 0.18-0.38 among the four groups. The K m values are three and four orders of magnitude higher in donors and contacts compared to males URT and males aerosol, respectively, which is due to the higher RNA values in the Sia et al. study [25].
For each group, we investigated the changes in viral RNA to infectious virus ratio, R u /V u , and found it to be time-dependent, which suggests that the two measurements are explaining different biological processes. Specifically, we found that the average R u /V u grows between 2 × 10 4 and 10 6 in the first 4 days in donors, between 2 × 10 4 and 3 × 10 8 in the first 5 days in contacts, between 8 and 600 in the first 3 days in males URT, and between 6 and 136 in the first 3 days in males aerosols, consistent with the data (Fig 7 blue, red, gold, and pink lines). Later R u /V u ratio predictions from models Eqs (3) and ( 4) are no longer reliable, as the infectious virus decays below limits of detection making the ratio unrealistically high.
In this study, we developed within-host and aerosol mathematical models of SARS-CoV-2 infection in golden hamsters that describe the kinetics of infectious virus and total viral RNA in the upper respiratory tract and aerosols. We fitted the models to data from two studies [25] and [24] which included four separate groups: three that were challenges with high viral dose (donors, males and females) and one that was challenged through infection following close proximity with an infected host (contacts). We estimated several key parameter values for each group and determined inter group variability. We found that the within-host basic reproductive number R 0 is less than one in all female hamsters, indicating limited viral spread. By contrast, R 0 in above one in all male hamsters from all groups (donors, contacts, males), indicating successful viral spread. Sex-based differences have been reported in other studies, which report that favorable outcomes in females are immune mediated [28,29]. The hamsters in the contact group had the lowest R 0 , averaging at 41, higher than other studies [13]. The other two groups, who were challenged with high viral dose, had slightly larger R 0 estimates, with average R 0 of 44 in donors and 60 in males. The death rate of productively infected cells ranged between 1.5 and 5.1 per day, corresponding to infected cells life-spans of 4 to 16 hours, similar within the groups and longer than in humans studies [13].
To model infectious virus shedding into the environment, we extended the upper respiratory tract model to account for infectious virus titer and total viral RNA emitted into aerosols. Within-host and aerosol model fitting showed that infectious virus titers get released into the air immediately after inoculation with peak shedding 21-37 hours after inoculation in males and 3-5 hours in females. Infectious viral shedding ends faster in females compared to males, with virus in aerosols losing replication-competency by 2.5-4.1 days post infection in females and 3.7-5.6 days post inoculation in males. The predicted loss of virus infectivity in aerosols before day six is consistent with experimental observations that have shown lack of transmission to naive hamsters co-housed with hamsters infected six days prior [25]. Our models and data herein predict that viral RNA persists in both upper respiratory tract and in aerosols long after replication-competent virus stops being detected, with RNA values staying above detectable levels at least a week after the infectious virus is lost. While in public health setting a SARS-CoV-2 diagnostic is determined by PCR assays, RNA levels are not always indicative of virus infectivity, with PCR specificity for detecting replication-competent virus decaying as the cycle threshold (Ct) values increase [30][31][32]32]. We used the models in this study to determine the connection between viral RNA and replication-competent virus levels, and found that the ratio of viral RNA to infectious viral titers is time-dependent and ranges between 10 2 and 10 8 RNA/TCID 50 in the first five days following infection, wider than in other studies [11,12]. We also found that after day five, the RNA to infectious virus ratio is no longer a reliable measurement of infectiousness, with the measured RNA values indicating the presence of genomic fragments, immune-complexed or neutralized virus, rather than replication-competent virus [33]. These estimates may be species and disease specific, with human studies reporting variable lengths of infectious viral shedding [16,32,34], with larger shedding windows during severe disease [34]. However, our results are consistent with the reported lack of transmission in contact golden hamsters co-housed with an infected donor at day six [25].
We have also investigated the population level relationship between the amount of RNA and the amount of infectious virus in a sample and found that the infectious virus increases in a density-dependent manner with the viral RNA, as suggested by previous work [13,24]. The results (consistent among the four groups considered) showed that when the viral RNA is high, the level of infectious virus saturates, consistent with our temporal results that show that it is unlikely we can predict infectiousness at high RNA:TCID 50 ratio. The turning point where viral infectiousness starts to saturate is study dependent, with differences due to experimental settings.
Our study has several limitations. We assumed that the initial target population is T(0) = 10 7 epithelial cells/ml, as in influenza [26]. This choice does not alter the results in this study, but it impacts the estimate of parameter p. We can obtain the same results for any other initial number of susceptible epithelial cells, as long as the p × T 0 is fixed at the current estimates. We assumed that both infectious virus clearance and RNA degradation rates are known, with infectious virus clearance set at influenza levels c + d = 10 per day (corresponding to life-span of 2.4 hours) and the degradation rate set arbitrarily at d = 1 per day (corresponding to lifespan of one day). Using sensitivity analysis we have found that changing the clearance rates to c + d = 15 and c + d = 5 does not influence the results (not shown). In all instances, however, the RNA degradation needs to be small, d = 1 or smaller, to explain the differences between infectious virus and viral RNA decay. Moreover, we had to include an additional removal of the aerosol infectious virus, which we assumed was due to infectious viral inactivation due to the elements, which we set at d 1 = d = 1 per day. We have also considered that all neutralized virus leads to RNA production. Further information is needed to determine the biological processes leading to increased degradation of infectious virus in aerosols compared to upper respiratory tract and those leading to the persistence of RNA in upper respiratory tract and aerosols after infectious virus is lost. Lastly, due to limited aerosol data in the females (with some subjects having just one data point above the limit of detection) we could not properly identify sex-specific differences and excluded this group from some of the analyses.
In conclusion, we have developed a within-host and aerosol model for SARS-CoV-2 infection in golden hamsters and used it to investigate the dynamics of viral RNA and infectious virus titers in URT and aerosols. We validated the models against data and used it to determine the temporal relationship between infectious virus, viral RNA and the probability of a nearby host getting infected. The results can guide interventions.
PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1009997 August 1, 2022