Introduction

Mathematical models have been extensively used to understand the epidemic characteristics of oubreaks, in predicting future outcomes, and in shaping the national responses regarding control measures [1, 2]. Despite the time pressure, a considerable amount of work has been dedicated to modeling the pandemic of novel coronavirus (SARS-CoV-2) infections that began in China in late 2019 [3–6]. Although most of these studies are based on existing epidemic models such as SIR and SEIR-models, several features of the COVID-19 pandemic have been independently explored, leading to different generalizations of similar dynamical models. On one hand, having a variety of models is central to get a notion of the model sensitivity, on the other, it shows that different assumptions are equally favorable to explain the observed data given the right set of parameter choices, whereas they might lead to different projections on how the epidemic would follow in the future [7, 8]. This variability in future projections becomes especially important when a perturbation, such as the imposition or release of the control measures, is introduced to the dynamical system.

A key epidemiologic variable that characterizes the spread of an infectious disease is the generation time [9], i.e., the time between successive cases in a chain of transmission. Li et al. [10] estimated the generation time distribution to have a mean of 7.5 (95%CI 5.5–19) days based on 6 observations, whereas Ganyani et al. estimated the generation time distribution to have a mean of 5.20 (95%CI 3.78–6.78) days for Singapore and 3.95 (95%CI 3.01–4.91) days for Tianjin [11], Bi et al. estimated the generation time distribution to have a mean of 6.3 (95%CI 5.2–7.6) days [12], He et al. estimated the generation time distribution to have a mean of 5.8 (95%CI 4.8–6.8) days [13], and Hiroshi et al. estimated the generation time distribution to have a mean of 4.7 (95%CI 3.70–6.00) days. Considering all these studies, infectiousness is estimated to decline quickly within 4 to 8 days on average.

Additionally, certain cases of transmission arouse concern about prolonged shedding of SARS-CoV-2 after recovery [14]. Moreover, several studies show proof of active virus replication in upper respiratory tract tissues and prolonged viral shedding even after seroconversion for COVID-19, implying that the contagious period of COVID-19 might last more than one week after clinical recovery in a fraction of patients [15, 16]. De Chang et al. reported patients to be virus positive even after the resolution of symptoms up to 8 days [17]. Similarly, Young et al. reported a median duration of 12 days for viral shedding [18], and Zhou et al. observed a median duration of 20 days [19]. Tan et al. reported a special case where the duration of viral shedding persisted for 49 days from illness onset [20]. Such examples indicate an uncertainty regarding the skewness of the generation time distribution. In additon to this uncertainty, several studies estimating the generation time suffer from short follow-up times, selection bias, and recall bias, which might miss the individual cases with prolonged shedding durations. Considering that the duration of infectiousness is a critical parameter in dynamical models used for predictive purposes, it is important to consider the epidemiological plausibility of a a more heavy tailed generation time distribution than the reported distributions in the literature and investigate its impact on model outcomes.

To do so, we first develop a generalized SEIR model by segregating the infectious compartment into two as “primarily infected” and “chronically infected” population. We assume that primarily infected individuals have a higher infectiousnesss within the time window conventially considered as the generation time, during when they have the potential to develop symptoms and therefore be hospitalized. Afterwards, we assume that the non-hospitalized infecteds transition to the chronically infected phase before recovery and become less infectious, but may stay infectious for a longer duration. By doing so, we include the possibility of a prolonged viral shedding window in our model. Individuals in the chronically infected phase are relevant both for diagnosis (a positive test result) and disease transmission, and we will explore the role of both aspects in explaining the observed data.

Using the incidence and fatality data from different regions of Italy and different states of the U.S., we first show that our model is also a plausible candidate for explaining the data observed prior to the easing of the lockdown measures (relaxation) for a variety of combinations of prolonged duration and level of infectiousness assumed for the chronically infected population. Based on this conclusion, we test the predictive power of different models using the daily confirmed cases data after the introduction of relaxation in Switzerland, including a model without the chronically infected population to represent the models conventionally used. Only Swiss data is used to test the predictive power of different models due to the public avaliability of different data types (estimates on the effective reproductive number, the number of daily confirmed cases, daily deaths, hospitalized and ICU patients) in high temporal resolution. Our results show that, in case of a gradual easing on the lockdown measures, the predictions of the model including the chronically infected population vary considerably from those obtained under a model in which prolonged infectiousness is not taken into account. This variability is especially important when national policies on control measures are being formed, and also for the healthcare systems if projections such as the occupancy of the hospital ward or the ICU are calculated using similar dynamical models.

Materials and methods

Mathematical model

To describe the dynamics of the COVID-19 pandemic, we generalize the susceptible-exposed-infected-removed (SEIR) compartmental model by including eight different states denoted by S(t), E(t), I_p(t), I_c(t), H(t), ICU(t), R(t), and X(t), representing the number of susceptible individuals, exposed (infected but not yet infectious) individuals, primarily infected individuals, chronically infected individuals, hospitalized patients, patients in ICUs, recovered (immune) individuals, and deceased individuals at time t, respectively. To model the prolonged viral shedding in case of COVID-19, we segregate the infectious compartment into two by introducing two different compartments, namely the primarily infected (I_p) and the chronically infected (I_c) individuals. After the incubation period is complete, exposed individuals become primarily infected where they stay infectious within the reported duration of the infectious period of COVID-19. Conventionally, these individuals are assumed to stop being infectious and therefore stop contributing to the disease spread when the generation time is complete. Our purpose by including another step before recovery, i.e., the chronically infected compartment, is to model a scenario such that the primarily infected individuals transition to a state where they are less infectious but they may stay infectious and be diagnosed for a longer duration than the generation time, i.e. continue spreading the infection with reduced transmissibility.

Transitions between different compartments are illustrated in Fig 1, which can be translated into a system of ordinary differential equations, where each arrow, i.e., each process, is associated with a rate. This system is given by the Eq set 1, including the rates of processes as model parameters, and describes the rate of change of compartments over time. Model parameters are given in Table 1 with their corresponding descriptions and prior distributions. An additional compartment C(t) is included in the Eq set 1 to calculate the cumulative number of the positively diagnosed cases in the community, and does not play any role in the disease dynamics. (1)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Illustration of the generalized SEIR model.

a) Notation of the compartments and their corresponding descriptions. b) Schematic of the dynamical model given by Eq set 1.

https://doi.org/10.1371/journal.pcbi.1008609.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Model parameters given with their descriptions, constrained ranges, and prior distributions.

https://doi.org/10.1371/journal.pcbi.1008609.t001

Time-dependent decrease in the transmission of SARS-CoV-2 due to lockdown measures is modeled by a sigmoid function [21], and denoted by r_lock(t), such that (2) where r_L, t_L, m_L, and s_L denote the final effect of the lockdown, start date of the lockdown, slope of the decrease in transmissibility, and the time delay between implementation and effect of the lockdown, respectively. r_lock(t) is used as a multiplicative factor in modeling the transmission rate in a time-dependent manner.

The reduced transmissibility of I_c is modeled via including a reduction coefficient r_c as a multiplicative factor to its transmission rate, representing the reduction in the infectiousness level of the primarily infected population when they move to the chronically infected phase. Introduction of r_c results in two different transmission rates β_p and β_c for I_p and I_c compartments, such that, (3) (4) (5) where , , 1/γ_p, and 1/γ_c denotes the basic reproduction number of the primarily infected population, the basic reproduction number of the chronically infected population, duration of primarily infected phase, and the duration of chronically infected phase, respectively. We assume that individuals who develop symptoms do so only during the primarily infected phase, and therefore hospitalization is only possible before they transition to the chronically infected phase. We do not assume any a priori information regarding the testing policy, therefore a positive diagnosis is possible for both primarily and chronically infected individuals, and they contribute to the cumulative number of the positively diagnosed cases with the rates and , respectively.

Results

Possibility of a chronically infected population

We find that the model that describes the data the best is dependent on the combination of the level and duration of the prolonged infectiousness, and the optimal choice of the {γ_c, r_c} combination varies among different regions and different countries. For Lombardy, as the duration of infectiousness becomes longer, models including the chronically infected population outperform the model without the I_c compartment (model w/o I_c) more often (Fig 2e), whereas the models with prolonged infectiousness (r_c < 100%) perform similarly to the model where the individuals can be diagnosed during the prolonged viral shedding window without being infectious (r_c = 100%) (Fig 2f). The absolute difference in RMSE values of the median of the model estimates differ by 30.7 daily confirmed cases (3.1% of the mean number of daily confirmed cases, S1 Table), and 6.17 daily deaths (3.4% of the mean number of daily deaths, S1 Table) the most when all r_c and γ_c values are considered (Fig 2c and 2d). The region of Emilia-Romagna presents a very similar behaviour to Lombardy (S1 Fig). In case of Piedmont and the state of New York, models with prolonged viral shedding (0% ≤ r_c ≤ 100%) outperform the model without the I_c compartment (model w/o I_c) more often, and models with lower levels of infectiousness (higher r_c values) clearly provide a better fit than the models with higher level of infectiousness (S2 and S5 Figs). For the state of Louisiana, all models perform similarly (S3 Fig). For the state of New Jersey, RMSE values for the number of daily confirmed cases are sensitive to the particular combinations of γ_c and r_c values (S4 Fig). Model without the I_c compartment (model w/o I_c) provides a better fit for very short and very long durations of prolonged infectiousness, and similar fits to the model with r_c = 100% for medium durations of prolonged infectiousness (S4 Fig). Maximum difference in the median RMSE values for the number of daily confirmed cases and the number of daily deaths for all combinations of r_c and γ_c values are provided in S1 Table, both in absolute values and relative to the mean of their corresponding data type. Parameter estimates with their corresponding means, standart deviations, and confidence intervals for all combinations of r_c and γ_c values are provided in the S2 Table.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Fitting and Root Mean Squared Error (RMSE) results for Lombardy.

Fitting and RMSE results for Lombardy, calculated using different levels and durations of infectiousness for the chronically infected population. Model outcomes (presented for 1/γ_c = 14 days) for the number of a) daily confimed cases and b) daily deaths using the data until the introduction of relaxation for model fitting, respectively. Darker shades of blue represent the fitting results with increased infectiousness of the chronically infected population, i.e., lower r_c values within the range 0 ≤ r_c < 100%. Fitting results for r_c = 100% are drawn in red, and the fitting results for the model without the I_c compartment (model w/o I_c) are drawn in pink. Data points that are used for fitting are drawn in black. Gray areas around the model outcomes represent the union of the 95% confidence intervals calculated for all models. RMSE values c) for the number of daily confirmed cases and d) the number of daily deaths for a given r_c and γ_c value used for fitting, where model w/o I_c represents the results for the model without the I_c compartment. e) Probability of the model without the I_c compartment (model w/o I_c) having a greater combined RMSE (CRMSE) value than the model with the I_c compartment for all levels of reduced infectiousness (r_c ≤ 100%) for different r_c and γ_c values. f) Probability of the model where individuals are being diagnosed without being infectious (r_c = 100%) having a greater combined RMSE (CRMSE) value than the model with individuals with a a prolonged infectiousness (r_c < 100%) for different r_c and γ_c values. Points in the gray areas represent the models that are providing a better fit more frequently than e) the model without the I_c compartment (model w/o I_c) and f) the model with r_c = 100%.

https://doi.org/10.1371/journal.pcbi.1008609.g002

Impact of relaxation

Data for Switzerland after the introduction of relaxation is used to test the predictive power of different models. As a demonstrative example, effect of the lockdown and the relaxation on infectiousness (r_lock(t) and r_relax(t)) are provided in Fig 3c and 3d for 1/γ_c = 14 days.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Fitting and prediction results for Switzerland.

a) Fitting results and relaxation predictions for Switzerland for the number of daily confirmed cases, calculated using different levels of infectiousness for the chronically infected population, assuming a duration of prolonged infectiousness of 1/γ_c = 14 days. Time dependent effects of the lockdown (r_lock(t)) and the relaxation (r_relax(t)) are illustrated in c) and d), respectively. Predictions drawn in darker shades of blue represent the fitting results with increased infectiousness of the chronically infected population, i.e., lower r_c values within the range 90 ≤ r_c < 100%. Fitting results for r_c = 100% are drawn in red, and the fitting results for the model without the I_c compartment (model w/o I_c) are drawn in pink. Data points that are used for fitting are drawn in black, and the data points used for comparing the predictive power of different models are drawn in green. e) RMSE values calculated using the prediction results for the number of daily confirmed cases for a given r_c and γ_c value, where model w/o I_c represents the results for the model without the I_c compartment. Models with the best predictive power (smallest RMSE value) are indicated by the bold black boxes. f) Probability of the model without the I_c compartment (model w/o I_c) having a greater RMSE value than the model with the I_c compartment for different levels of reduced infectiousness (90% ≤ r_c ≤ 100%) and γ_c values, calculated over the predicted data points. g) Probability of the model where individuals are being diagnosed without being infectious (r_c = 100%) having a greater RMSE value than the model with individuals with a a prolonged infectiousness (r_c < 100%) for different r_c and γ_c values, calculated over the predicted data points. Points in the gray areas represent the models that are providing a better fit more frequently than f) the model without the I_c compartment (model w/o I_c) and g) the model with r_c = 100%. b) Prediction results of the first 5 models (r_c = {93%, 94%, 95%, 96%, 97%}) with the lowest RMSE value (best predictive power), model with r_c = 100%, and model without the I_c compartment (model w/o I_c) for + T days into the future from the last data point observed, where 1/γ_c = 14 days. Predictions for the model with the best predictive power (r_c = 95%), for the model with r_c = 100%, and the model without the I_c compartment (model w/o I_c) are highlighted in blue, red, and pink, respectively.

https://doi.org/10.1371/journal.pcbi.1008609.g003

We observe that all models provide almost identical fits for the data prior to the introduction of relaxation (S6 Fig), but they substantially differ in their predictions regarding after relaxation even for small differences in the infectiousness levels (r_c values) for the chronically infected population (Fig 3a and 3b, demonstrative example for 1/γ_c = 14 days). As the predicted point moves further in time, the quantitative difference between the predictions of different models deviate from each other even more. For 60 days after the last data point observed, the model with the best predictive power (model with r_c = 95%) predicts a median of 1439 daily confirmed cases, whereas model with r_c = 100% predicts a median of 126, and model without the I_c compartment (model w/o I_c) predicts a median of 190 daily confirmed cases (Fig 3b, T = 60), indicating a discrepancy by one order of magnitude. Similar to the results provided for the first stage of fitting, we find that it is more advantageous to use a model including the chronically infected population with low levels of infectiousness (Fig 3f and 3g for estimates bootstrapped within 50%, and S7 Fig for estimates bootstrapped within 95% confidence intervals). The model providing the lowest RMSE value between the predictions and the data is always a model with prolonged infectiousness for a wide range of γ_c values (Fig 3e). However, note that these results are valid under our particular choice of parametrization of the relaxation dynamics, and the resulting relative change in the transmissibility during the relaxation period.

The fact that observed data can be explained equally well by various combinations of r_c and γ_c values is partially due to the flexibility of the fitting procedure, which allows other parameters to be adjusted for a given {r_c, γ_c} pair. Most parameters are free to vary, but their prior distributions are informed such that the hyperparameters (parameters of the prior distributions) align with the reported values in the literature (Table 1). As an example, both the incubation period (1/τ) and the duration of infectiousness of the primarily infected population (1/γ_p) have the mean of 2.5 days, resulting in a generation time distribution with a mean of 5 days, in agreement with the reported values in the literature for COVID-19 (see Introduction). Similarly, the basic reproduction number of the primarily infected population is normally distributed with a mean of 2.5, which is the average value reported for basic reproduction number of COVID-19 in many countries [10, 27]. Mean values of the prior distributions of the parameters related to hospitalization (γ_H, γ_ICU, ϵ_H, ϵ_H2I, and ϵ_x) are adopted from Ferguson et al. [28] and Verity et al. [29], and given a variance such that they can be adjusted specifically for each country during the fitting procedure.

The median of the posterior distributions for , R₀, and r_L provide a good example to demonstrate the flexibility of the fitting procedure (Fig 4). As expected, the basic reproduction number of the primarily infected population () (Fig 4a), the total basic reproduction number () (Fig 4b), and the final reduction in infectiousness due to lockdown (1 − r_L) (Fig 4c) are estimated to be lower for a given duration of infectiousness (1/γ_c) as the infectiousness of the chronically infected population decreases (as r_c increases) to explain the observed data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Median of the posterior distributions of parameters.

Median values of the posterior distributions of the a) basic reproduction number of the primarily infected population (), b) total basic reproduction number (), and c) the final reduction in infectiousness due to lockdown (1 − r_L), for a given {r_c, γ_c} pair.

https://doi.org/10.1371/journal.pcbi.1008609.g004

Parameter estimates for the Swiss data with their corresponding means, standart deviations, and confidence intervals for all combinations of r_c and γ_c values are provided in the S2 Table.

Discussion

The model presented in this work explores the plausibility of an epidemiological model with a prolonged viral shedding window for the COVID-19 pandemic, and investigates both its impact and predictive capabilities on the outcomes of a gradual easing on the lockdown measures (relaxation) given different assumptions on the infectiousness level and duration of a chronically infected population.

Our results show that including a chronically infected population, i.e., individuals that are less infectious but infectious for a longer duration, is not a possibility that can be easily rejected from an epidemiological perspective. This conclusion is based on two main results. First, neither the presence nor the absence of chronic transmission is identifiable from population-level data. The data that has been observed until relaxation can be explained equally well by the model with prolonged viral shedding for a variety of different levels and durations of prolonged infectiousness as by the model without prolonged viral shedding. Although this is partially due to the flexibility of the fitting procedure, the choice of hyperparameters (parameters of the prior distributions) indicates that all fits for a given infectiousness value are possible for a set of reasonable model parameters, and therefore as favorable as the conventional models from a modeling perspective.

Second, even if the presence of a chronically infected population cannot be proven, its introduction to the model structure has a considerable impact on the relaxation outcomes. In case of a gradual easing on the lockdown measures, the predictions of the model including the chronically infected population vary considerably from those obtained under a model in which prolonged infectiousness is not taken into account. Although the level of infectiousness might be low, its impact during the prolonged viral shedding window is significant in terms of predicting the outcomes of a gradual relaxation, indicating that even small differences in prolonged infectiousness levels might change the course of an epidemic when they are present for a certain duration. This is especially important for the healthcare systems if projections such as the occupancy of the hospital ward or the ICU are calculated using similar dynamical models.

The fact that observed data can also be explained with a model including prolonged viral shedding raises certain questions about the interpretation of the epidemic curve, acquired immunity, and the current testing policies. Assuming a relatively short generation time for a model that does not consider a prolonged viral shedding window results in more optimistic projections about epidemic control, as clearly demonstrated in Fig 3. Based on this, countries that were very successful in their initial control measures and therefore experienced a very steep decline in the number of daily confirmed cases might choose to ease the control measures too soon. We still lack a full understanding of the viral shedding window of COVID-19, and therefore might have a biased opinion on the number of infectious individuals in the community. This once again emphasizes the infectiousness of COVID-19 and the significance of frequent testing although the number of cases are in decline.

Using simplified compartmental models such as the one in this study has certain limitations. First, it does not consider the stochastic effects that the system is subject to, which become more important as the number of infecteds decrease in the community. Second, it assumes a well-mixed population, and does not consider the contact structure and the demographic information which are both relevant to the disease spread. Nevertheless, we believe that these two drawbacks of our modeling approach influence all models with and without the prolonged viral shedding to a similar degree, if not penalizing the models with prolonged viral shedding for producing more pessimistic projections since the number of infecteds will be higher in frequency relative to the model without the chronically infected compartment. Additionally, the standard SEIR model assumes constant rates of transition between the exposed, infectious, and recovered classes, leading to waiting times that individuals spend in these states being exponentially distributed [30], resulting in an exponential distribution on the generation time as well. Although mathematically convenient, this assumption is shown to be epidemiologically unrealistic, and less dispersed distributions such as gamma distribution should be used instead [31, 32].

Compartmental model structures are based on the underlying epidemiological and demographic interactions of a particular disease. Given that there are many choices for these interactions, the number of possible combinations are enormous [33, 34]. Our choice of including a chronically infected comparment in the model structure was inspired by the evidence indicating an uncertainity regarding the generation time distributions, but this approach is only one way to extend the basic SEIR model for the COVID-19 pandemic. There are still several open questions regarding the transmission dynamics of COVID-19, meaning that there are many other alternatives of modification a modeler could consider depending on the research question in hand. These alternatives are also potential candidates which would describe the data equally well, and offer reasonable predictions.

One methodological limitation of pure model fitting is the parameter identification problem, especially in the early stages of an epidemic [35]. As clearly demonstrated in Fig 4, models with different assumptions on the duration and the level of prolonged infectiousness lead to equally good descriptions of the observed data by adjusting the parameter values accordingly. Therefore, even if a prolonged viral shedding window exists for COVID-19, it would not be possible to quantify the precise level or the duration of prolonged infectiousness by using model fitting purely. Ultimately, these quantities should be measured or estimated from the relevant type of data.

Another potential limitation is the dependency of the goodness of fit of a given model on the quantification of the impact of relaxation, which inevitably affects the model selection procedure. Although we do not perform any fitting on the data belonging to the relaxation phase, we indirectly inform our predictions by shaping the change in transmissibility (β) via using the normalized values of the R_e estimates provided by the Swiss National COVID-19 Science Task Force, which in turn are calculated using the data on the number of daily cases. Different assumptions on the change of transmissibility during relaxation might alter the infectiousness level that is optimal in predicting the relaxation outcomes. With that being said, our results suggest that the model without the chronically infected compartment heavily underpredicts the case numbers, as clearly seen in Fig 3. Although it is still debated whether the patients who recover from COVID-19 and test positive for the virus after their recovery are still infectious or not, it is clear that these positive test results contribute to the data on the number of daily confirmed cases. However, current modeling studies regarding COVID-19 neglect this fact and assume that the probability of detecting an infection decreases strongly after the mean generation time. Our results show that this assumption might lead to an underestimation of both the reproduction number and the effect of the lockdown (Fig 4), leading to a potential underprediction for the relaxation outcomes.

In conclusion, it is not possible to either prove or disprove the existence of a compartment of chronically infected individuals purely by modeling based on epidemiological data. Our results only provide tentative evidence to consider a chronically infected population as an alternative modeling approach in addressing the knowledge gap on the transmission dynamics of COVID-19. Such an hypothesis must be tested by incorporating data regarding the timing of transmission events, contact histories, and corresponding test results. Furthermore, more clinical and virological diagnostic studies are necessary to establish the biological links between viral load, active viral replication at different sites of the body, severity of symptoms, and a positive test result to infer the infectiousness of an individual over time. Including a chronically infected population in our model was motivated by the evidence reported for prolonged viral shedding in the literature [14–20], and attempted to test whether this is also a plausible descriptive and predictive modeling approach. Given that different assumptions on the infectiousness duration and level during a prolonged viral shedding window can result in similar descriptions of the observed data prior to the introduction of relaxation, and large differences of epidemic projections after relaxation, it is important to consider a chronically infected population from a modeling perspective when national policies are being imposed.

Supporting information

Acknowledgments

We gratefully acknowledge Dr. Julien Riou for his valuable comments and discussions.