Isolation may select for earlier and higher peak viral load but shorter duration in SARS-CoV-2 evolution – Nature Communications

Characterizing the time-series pattern of viral load for SARS-CoV-2 variants

We analyzed the viral load of SARS-CoV-2 for the pre-Alpha (non-VOI/VOCs)^14,31, Alpha^31,32, and Delta variants^31,32 using a simple mathematical model describing viral dynamics^15,16,17,18 (see Table S1 and METHODS). To consider inter-individual variation in the patients’ viral loads, we employed nonlinear mixed effect models to estimate the parameters. The estimated population parameters are shown in Table S2. The typical behavior of the mathematical model, Eqs. (4, 5), using these parameters is shown in Fig. 1a and Fig. S1 for the pre-Alpha (black), Alpha (blue), and Delta variants (red), respectively. Note that the time of infection (i.e., t=0) is estimated in our model fitting by using fixed initial values³³ (see METHODS). These quantified virus infection dynamics highlight how the properties of the virus differ variant by variant.

To characterize and compare the infection dynamics of these SARS-CoV-2 variants, we quantified the distributions of the peak viral load, peak time, and duration of viral shedding, D, (Fig. 1b–d). Here duration of viral shedding indicates the time for which the viral load is detectable. We found that peak time did not differ significantly between the pre-Alpha and Alpha variants (Fig. 1c). The peak viral load also was not significantly different between the Alpha and Delta variants (Fig. 1b, d). However, interestingly, we found that the peak viral loads for the Alpha and Delta variants were higher than for the pre-Alpha variants (\(p=2.9\times {10}^{-7}\) and \(8.8\times {10}^{-6}\) by the Mann-Whitney-Wilcoxon test, respectively) (Table 1). The peak time and duration of viral shedding for the Delta variant were shorter than those for the pre-Alpha and Alpha variants (for the peak time: p-values are less than \(2.2\times {10}^{-16}\) for the duration of viral shedding: \(p=1.4\times {10}^{-15}\) and \(1.1\times {10}^{-6}\) by the Mann-Whitney-Wilcoxon test, respectively) (Table 1). Taken together, our results highlight that the duration of viral shedding was shortened and the peak viral load was increased and advanced as SARS-CoV-2 evolved. This implies that SARS-CoV-2 is evolving to a more “acute phenotype”³⁴, mainly characterized by an advanced peak time and a shorter duration of viral shedding compared with the pre-Alpha variants.

To further characterize the time-series pattern of viral load for SARS-CoV-2 variants during evolution of the virus (i.e., the Alpha variant emerged during the spread of the pre-Alpha variants, and the Delta variant emerged after the spread of the Alpha variant), we first investigated the relation between the duration of viral shedding, D, and the cumulative log-transformed viral load, \({V}_{{{{{{\rm{total}}}}}}}\), defined as the area under the curve (i.e., AUC [\({\log }_{10}{{{{{\rm{viral}}}}}}\; {{{{{\rm{RNA}}}}}}\; {{{{{\rm{copies}}}}}}/{{{{{\rm{ml}}}}}}\times {{{{{\rm{days}}}}}}\)]) (Fig. 1e). We found \({V}_{{{{{{\rm{total}}}}}}}\) increased as D increased (i.e., positive correlation) (see METHODS). Interestingly, there was also a trend that \({V}_{{{{{{\rm{total}}}}}}}\) decreased as SARS-CoV-2 evolved: \({V}_{{{{{{\rm{total}}}}}}}\) for the Delta variant was smaller than that for the pre-Alpha and Alpha variants (\(p=1.1\times {10}^{-7}\) and \(9.4\times {10}^{-6}\) by the Mann-Whitney-Wilcoxon test, respectively). In addition, to understand the epidemiological consequences of these variants, we also calculated the transmission probability of each variant depending on the time after infection (see Fig. 1f and METHODS). We found that the transmission probability of the Alpha variant remained a high transmission probability for a longer period of time than the pre-Alpha variant. Furthermore, the transmission probability of the Delta variant peaked at 3.0 days after infection⁹ but decreased to almost 0 by 8.0 days, suggesting tFhe evolutionary process promotes effective and rapid transmission within human populations (see below for details).

Prediction on time-series pattern of viral load of SARS-CoV-2 evolution with NPIs

In the comparison of viral shedding dynamics, the Alpha variant showed a higher peak viral load than the pre-Alpha variants (Fig. 1b). Although the duration of viral shedding for the Alpha variants was less than that for the pre-Alpha variant (Fig. 1d), the transmission probability of the Alpha variants remained higher for a long time (Fig. 1f). Note that a similar difference between the pre-Alpha and Delta variants can be observed (Fig. 1b, d, f). These findings demonstrate that the Alpha-like and the Delta-like variants showed a high peak viral load but a relatively short duration of viral shedding as the virus evolved among human populations.

In contrast, in a comparison of the Alpha and Delta variants, the selective force to facilitate the evolution of the Delta-like variant is not trivial. This is because the peak viral load and the duration of viral shedding were similar between the Alpha and Delta variants (Fig. 1b, d), and therefore the highest transmission probability maintains almost the same level and length for these variants (Fig. 1f). The difference between the Alpha and Delta variants is the peak time, that is, the peak time of the Delta variant was significantly earlier than for the Alpha variant (Fig. 1c). Our hypothesis is that the Delta-like variants showing an earlier peak time but similar peak viral load and duration of viral shedding (i.e., acute phenotype) had a general advantage under human-mediated selection pressures (i.e., earlier infections result in more descendants). COVID-19 patients show both symptomatic and asymptomatic infection, and the incubation period and the proportion of asymptomatic patients are decreasing and increasing, respectively, as SARS-CoV-2 evolves^35,36. To prevent the chain of transmission, patients showing symptoms are usually isolated or hospitalized or subjected to other NPIs^37,38. Symptomatic patients primarily generate secondary infections prior to symptom onset, that is, during the pre-symptomatic period. On the other hand, asymptomatic patients can be infectious during the whole duration of viral shedding. Below we show that a strong NPI, defined as the isolation of symptomatic infected individuals after symptom onset, could be a driving force for SARS-CoV-2 evolution. Note that although we can easily incorporate variant-specific virological characteristics such as infectivity per virus by changing the probability of transmission (see METHODS), we ignore these differences here for simplicity.

NPI-driven SARS-CoV-2 evolution

Using the probabilistic multi-level model, which considers both individual-level virus infection dynamics and population-scale transmission, as shown in Eqs. 1–3, 6, and the genetic algorithm (GA), which mimics virus phenotypic evolution³⁴, we explored how the time-series patterns of viral load evolved under NPIs. The schematic diagram of our multi-level modeling is described in Fig. 2, and a detailed explanation is provided in the METHODS.

A schematic of the multi-level population dynamics model is depicted. Virus transmission occurs from an infected individual to susceptible individuals depending on the infected individual’s infectivity, which depends on their viral load. At each time step (i.e., day), the focal infected individual has contact with multiple susceptible individuals. Here the contact numbers are assumed to follow a negative binomial distribution, which does not depend on the viral load. The shape parameters k and θ for this binomial distribution are presented in Table S3. The sum of newly infected individuals (i.e., [\(c\left(t\right)\): number of contacted individuals per day (everybody is susceptible)] × [\(\rho \left(t\right)\): transmission probability per contacted individual]) during the infectious period is calculated as \({R}_{{TP}}\), called the “transmission potential”. Note that the focal infectious individuals are assigned one of two properties: symptomatic or asymptomatic. Whether being symptomatic or asymptomatic is randomly assigned by a constant probability f or 1-f, respectively. It is assumed that the transmission chains from the focal symptomatic infectious individuals are perfectly inhibited by isolation after symptom onset, \({T}^{*} < t\). Otherwise, \({R}_{{TP}}\) is calculated using the whole duration of viral shedding of the focal asymptomatic infectious individual. Thus, \({R}_{{TP}}\) for the symptomatic and asymptomatic individuals are different here.

Briefly, from the assumption that the virus evolves to increase its “transmission potential”, defined as \({R}_{{TP}}\) (i.e., the total number of secondary cases generated throughout the infectious period), we calculate \({R}_{{TP}}\) under various patterns of infection rate, β, and virus production rate, p. Thus, the transmissibility fitness landscape of \({R}_{{TP}}\) is constructed as a function of \(\left(\beta,p\right)\), given the incubation period (i.e., time from infection to symptom onset), \({T}^{*}\), and the proportion of symptomatic patients, f. Note that the transmissibility fitness landscape of \({R}_{{TP}}\) is constructed as a function of only β and p, because the other parameters or initial values are estimated or determined (see METHODS). Under this analysis, we considered the (additional) effect of isolation, which is only applicable to symptomatic individuals and which assumes that isolation always perfectly prevents the transmission chain (i.e., perfect isolation; \({T}^{*} < \, t\)). In other words, the (baseline) effect of isolation, which is applicable to all individuals regardless of symptomatic or asymptomatic infection (e.g., wearing a face mask, social distancing), is assumed to be involved in the parameters in our multi-level model without loss of generality.

In the absence of symptomatic infection (i.e., all patients are asymptomatic: f = 0), we analyze SARS-CoV-2 evolutionary dynamics by exploring the optimal set of β and p, which characterizes the time-series pattern of viral load, with GA (see Table S3, S4 and Algorithm S1). Without symptomatic infection, there is no effect of symptom-dependent isolation on the transmission chain at all. The optimal sets of \(\left(\beta,p\right)\), the transmissibility fitness landscapes, the trajectories of \({R}_{{TP}}\) along the course of GA, and its distribution are described in Fig. S2a–c, respectively. The time-series patterns of viral load with the optimal parameters and the corresponding timing of peak viral load (i.e., peak time) are described in Fig. S2d, e, respectively.

In contrast, in the absence of asymptomatic infection (i.e., all patients are symptomatic; the fraction of symptomatic infection is \(f=1\)), all individuals lose their transmissibility by NPIs after symptom onset (\({T}^{*} < t\)). Depending on the incubation period, \({T}^{*}\), the distributions of optimal sets of \(\left(\beta,p\right)\) are located in different parameter ranges, and the corresponding transmissibility fitness landscapes are altered (Fig. S3a, b). Because isolation can prevent the transmission chain before the focal infected individuals became highly infectious, given a small \({T}^{*}\) such as \({{T}}^{{*}}{=}{1}{.}{0}\) day (i.e., almost no transmissions occur during the pre-symptomatic period), \({{R}}_{{{TP}}}\) with the optimal parameters is basically less than 1 (Fig. S3b). On the other hand, given a large \({{T}}^{{*}}\) (e.g., \({{T}}^{{*}}{=}{10}{.}{0}\) days), isolation is delayed, and therefore the transmissions are allowed during the pre-symptomatic period. Thus, \({{R}}_{{{TP}}}\) increases as \({{T}}^{{*}}\) increases (Fig. S3b). Note that the effect of isolation is weak for \({{T}}^{{*}}{=}{10}\) because the viral load already drops to a very low level at the moment when isolation is applied to the symptomatic patients, and thus the transmission probability is low. In fact, the virus evolution shows similar patterns observed in the absence of symptomatic infection (\({f}{=}{0}\)) (Fig. S2a–c vs. Fig. S3a–c for \({{T}}^{{*}}{=}{10}{.}{0}{,}\) for example). The time-series patterns of viral load and the corresponding timing of peak viral load are also described in Fig. S3d, e, respectively. Interestingly, compared with the case where all patients are asymptomatic (f = 0), acute phenotypes showing an earlier peak but similar peak viral load and duration of viral shedding are selected (i.e., isolation-driven evolution). Thus, in the absence of asymptomatic infection, selection generated by isolation is maximized and drives convergent evolution of viral load dynamics (Fig. S3d). Additionally, because of smaller transmission potential (i.e., secondary [focal] infectious individuals may not be effectively generated and sometime less than one), isolation-driven evolution might be unlikely to occur under this situation.

In the realistic situation, there are both asymptomatic and symptomatic infections (i.e., \({0}{ < }{f}{ < }{1}\)), and only symptomatic patients interfere with secondary infection by isolation after symptom onset. We next evaluated virus evolution under the different incubation period of \({{T}}^{{*}}\) with the fixed fraction of \({f}{=}{0}{.}{3}\) (i.e., 30% of infected individuals show symptoms). Note that there is substantial uncertainty in the value of \({f}\), and we conduct supplementary analyses for different values of \({f}\) (see below). As described in Fig. 3a, the optimal sets of \(\left({\beta }{,}{p}\right)\) are widely distributed for \({{T}}^{{*}}{=}{1}{.}{0}\) and \({3}{.}{0}\) days, but those are narrowed for \({{T}}^{{*}}{=}{6}{.}{0}\) and 10.0 days. This is because the transmissibility fitness landscapes of viruses differ between asymptomatic and symptomatic patients in the case of relatively small incubation periods, whereas the difference is small for large periods (Fig. 3b). We note that having a higher peak viral load while maintaining a rather long duration favors an increase in \({R}_{{TP}}\), and viral dynamics with this characteristic are achieved when β is low and p is relatively high, such as optimal sets of \(\left(\beta,p\right)\) for the asymptomatic case (or in the absence of isolation) in Fig. 3b.

a Genetic algorithm (GA) exploring the evolutionary trajectories on the (β,p) plane until the generation of 300 is applied, depending on different values of the incubation period, \({T}^{*}\). Here β andp are the infection rate and the virus production rate, respectively. All symptomatic individuals lose their transmissibility by isolation after symptom onset (\({T}^{*} < t\)). The white dots represent the endpoint of 100 independent simulation runs, and the contour lines are the kernel density estimation of their distribution. The colored dot in each panel is the mean value of the white dots, which represents the set of evolutionary outcomes of (β,p) under the parameters we used. The colors used here are independent to them in Fig. 1. The black line is the mean trajectory of the GA through 300 generations. b The mean transmissibility fitness landscapes aggregated solely from the asymptomatic (top row) and symptomatic (bottom row) individuals are described, respectively, using 100 runs of GA. The white dot represents the maximum value of the mean transmissibility fitness, \({R}_{{TP}}\). Note the asymptomatic individuals have larger transmissibility fitness than symptomatic cases on average, regardless of \({T}^{*}\) values. c The trajectories of \({R}_{{TP}}\) along the course of GA with different \({T}^{*}\) are calculated. The gray dotted lines are the mean trajectory over 100 trials of colored lines. d The time-series patterns of viral load with the optimal parameters of (β,p) with different \({T}^{*}\), which were obtained in a, are shown. e The contour-plot for the timing of peak viral load (i.e., peak time) is shown. Each curve is colored accordingly. The gray region is the parameter range satisfying \({R}_{{TP}} \, < \, 1\).

In the search for the optimal parameter set by GA, depending on the type of focal infectious individual (who is stochastically chosen), directions of virus evolution stochastically change so that neither peak of the transmissibility fitness landscape can be unreachable (compare the colored dots in Fig. 3a with the white dots in Fig. 3b for \({{T}}^{{*}}{=}{1}{.}{0}\) and 3.0) unless the evolutionary rate is extremely fast. Thus, an intermediate parameter set becomes the equilibrium state. Moreover, we calculated \({{R}}_{{{TP}}}\) for each infected individual: the trajectories of \({{R}}_{{{TP}}}\) with the optimal parameters vary around 1, especially with short incubation periods (Fig. 3c), because \({{R}}_{{{TP}}}\) depends on the types of focal infectious individual and \({{R}}_{{{TP}}}\) for the symptomatic individuals is basically smaller than that for the asymptomatic individuals (sometimes less than 1) because of isolation (Fig. 3b).

As observed in the absence of asymptomatic infection, the time-series patterns of viral load (Fig. 3d) with the optimal parameters (Fig. 3e) also depend largely on the incubation period \({{T}}^{{*}}\) because of isolation. Of note, in the realistic situation with \({0}{ < }{f}{ < }{1}\), because (secondary) infectious individuals are mainly maintained via transmissions by asymptomatic individuals, isolation-driven evolution for relatively small incubation periods may occur, and the peak viral load is advanced. Fig. 3d, e also demonstrate that faster peak time is associated with larger value of the product of p and β. This is due to the combined effect of larger p and β values, resulting in both rapid and increased amount of virus production during an early phase of infection. For our sensitivity analysis, we conducted the same simulations for \(f={{{{\mathrm{0.5,0.6,0.7}}}}}\) and obtained similar conclusions (see Fig. S4).

We have thus far assumed ‘perfect isolation’ in our results, where transmission is completely halted once patients become symptomatic. However, we understand that in practice, isolation is sometimes imperfect, leading to continued virus transmission. To account for this, we introduced a factor ξ, modifying the contact numbers post symptom onset (\({C}^{*}\left(t\right)=C\left(t\right)\times \xi\)), with ξ representing incomplete isolation probabilities of 0.1 and 0.3. In general, this extension yielded consistent outcomes with our initial findings (Fig. S5). However, it indicated that incomplete isolation leads to a less pronounced evolution towards an earlier, higher peak in viral load dynamics at a small \({T}^{*}\) (such as \({T}^{*}=1\): it is unreasonably small in the context of COVID-19 and therefore we may exclude this situation thought) (Figs. S5–1). Even in cases where \({T}^{*}\) is small, and it significantly drives virus evolution, the deterministic selective pressure is weakened due to the probabilistic failure of isolation. This situation is analogous to a scenario where all patients have asymptomatic infections as a result, making advancing the peak viral load no longer beneficial because of the short pre-symptomatic period. The trend becomes more prominent when the probabilistic failure of isolation occurs frequently (e.g., \(\xi=0.3\) in Figs. S5–2).

As we confirmed in Figs. 3, Fig. S2, S3, and S4, the higher proportion of symptomatic infection enhances the selection pressure for virus evolution via isolation. Taken together, these results suggest that while isolation inhibits the chain of transmission from symptomatic focal infectious individuals after symptom onset, in the realistic situation, SARS-CoV-2 evolves toward the acute phenotypes.

Effect of prior immunity on NPI-driven SARS-CoV-2 evolution

Considering another modeling assumption, our initial mathematical model did not account for the influence of COVID-19 vaccinations or prior infections. Given the widespread infection and vaccination leading to high-level population immunity worldwide, this factor could have a profound impact on the evolution of the virus. Some studies suggest that COVID-19 vaccinations may affect the SARS-CoV-2 infection dynamics^32,39. Accordingly, we conducted further analysis to investigate the impact of prior immunity caused by vaccinations or prior infections on our findings.

We first reanalyzed SARS-CoV-2 infection dynamics in Fig. S6, including additional 68 participants (14 and 54 patients for the Alpha and the Delta variants, respectively), who were excluded due to prior infections, and obtained the consistent statistical trends in line with our findings described in Fig. 1a, reinforcing the robustness of our conclusions. Specifically, we observed that the duration of viral shedding keeps decreased in evolution. The peak time did not exhibit significant differences between the pre-Alpha and Alpha variants, and the peak viral load did not significantly differ between the Alpha and Delta variants. However, the peak viral loads for both the Alpha and Delta variants were higher compared to the pre-Alpha variants (\(p=0.002\) and \(3.3\times {10}^{-7}\), respectively). Moreover, Delta variant exhibited a shorter peak time in comparison to both the pre-Alpha and Alpha variants (\(p=2.2\times {10}^{-16}\) for both cases) (Fig. S6). These confirmed that the statistical trend for SARS-CoV-2 evolution is remained regardless of the prior infections.

Furthermore, we extended our in silico analysis to incorporate a vaccinated population, considering vaccination rates set probabilistically at 50% and 80%. For those vaccinated, we increased the clearance rate (δ) by 1.5 times. This assumption for clearance rate shortened the viral shedding period by about 3 days, and this is larger than the previously reported shortening period of 2 days³². Thus, a 1.5-fold increase in clearance rate is sufficient to account for an increase in clearance rate caused by the vaccination. The incorporation of this stochastic environment led to a modest increase in the variance of optimal virus parameters. However, it is essential to note that these changes in the model maintained our main conclusion. The mean values derived remained consistent with our original findings: viruses were selected for had an earlier peak and higher viral load dynamics, but a shorter duration of infection (Fig. S7). The selective pressure for increased transmissibility shapes the viral load dynamics, and isolation measures are likely to be a critical driver of these evolutionary transitions even considering the prior immunity.

Validating SARS-CoV-2 evolution with Omicron variants

The Omicron variant (originally B.1.1.529) was first isolated in South Africa in October 2021^10,40. After that time, of the identified Omicron subvariants, the BA.1 subvariant rapidly spread worldwide and became the most prevalent SARS-CoV-2 variant in many countries (Fig. 4a). Here, we consider simply the situation that BA.1 emerged during the spread of the Delta variant, which was the most predominant variant worldwide before BA.1, under NPIs.

a An overview of frequency changes in SARS-CoV-2 variants based on data from CoVariants.org is shown. The left and right panels show the trend in the United States and the United Kingdom, respectively, as examples (because the data we used here were from the USA and UK). See the labels for each color curve corresponding to each VOC. b The estimated viral dynamics using our model and publicly available data for the BA.1 subvariant (green) are described along with the Delta variant (red) using the best-fit population parameters. The solid curves correspond to the solution of Eqs. (4, 5) using the best-fit population parameters, and the error bands represented as shadow regions indicate the 95% interquantile range of the inferred viral load interval for all individuals at each time point. c The distributions of the peak viral load, peak time, and duration of viral shedding for each patient are shown with their mean values (vertical solid bars). A statistically significant difference between two groups was found only for peak viral load (p-values by the two-sided Mann-Whitney-Wilcoxon test are \(1.3\times {10}^{-6}\) for peak viral load, 0.09 for peak time, and 0.48 for duration).

To validate our hypothesis that isolation may be a selection pressure to further drive the evolution of viral phenotypes, we additionally used the longitudinal viral load of SARS-CoV-2 from the BA.1 subvariant in 49 infected patients⁴¹ (Table S1), and similarly analyzed those data together with the same data for the Delta variant (see METHODS in detail). The individual viral loads for the Delta variant (red) and BA.1 subvariant (green) are described in Fig. S8, and the estimated population parameters are shown in Table S5. We found that the time-series pattern of viral load for the BA.1 subvariant differed from that for the Delta variant as described in Fig. 4b. We further characterized and compared the distributions of the peak viral load, peak time, and duration of viral shedding in Fig. 4c. Interestingly, while we found almost the same peak time for the BA.1 subvariant viral load with a similar duration of viral shedding (p = 0.09 and 0.48 by the Mann-Whitney-Wilcoxon test, respectively), the comparison of peak viral load suggested a lower peak for the BA.1 subvariant than for the Delta variant (\(p=1.3\times {10}^{-6}\) by the Mann-Whitney-Wilcoxon test). On the other hand, as we expected, the BA.1 subvariant showed an earlier peak time and shorter viral shedding compared with the pre-Alpha and Alpha variants (see Fig. 1a–d vs. Fig. 4b, c). Altogether, these results imply that SARS-CoV-2 evolution toward an earlier peak viral load was maintained but slowed down. Interestingly, similar evolutionary trends were confirmed in a non-human primate model infected with different VOCs of SARS-CoV-2⁴², obviously without NPIs in these cases. We discuss a mechanism behind those convergent phenotypic changes in detail later (see Discussion). Our quantitative empirical clinical data analysis based on our hypothesis could support that SARS-CoV-2 variants have evolved their acute phenotypes via a human-mediated selection pressure such as isolation.