Predictability of viral load kinetics in the early phases of SARS-CoV-2 through a model-based approach

Abstract

A pipeline to evaluate the evolution of viral dynamics based on a new model-driven approach has been developed in the present study. The proposed methods exploit real data and the multiscale structure of the infection dynamics to provide robust predictions of the epidemic dynamics. We focus on viral load kinetics whose dynamical features are typically available in the symptomatic stage of the infection. Hence, the epidemiological evolution is obtained by relying on a compartmental approach characterized by a varying infection rate to estimate early-stage viral load dynamics, of which few data are available. We test the proposed approach with real data of SARS-CoV-2 viral load kinetics collected from patients living in an Italian province. The considered database refers to early-phase infections, whose viral load kinetics are not affected by mass vaccination policies in Italy. Our contribution is devoted to provide an effective computational pipeline to evaluate in real time the evolution of infectivity. Comprehending the factors influencing the in-host viral dynamics represents a fundamental tool to provide robust public health strategies. This pilot study could be implemented in further investigations involving other respiratory viruses, to better clarify the process of viral dynamics as a preparatory action for future pandemics.

Keywords: Mathematical epidemiology; SARS-CoV-2 viral load; SIR with age of infection;
                  Uncertainty quantification methods.

AMS Subject Classification: 92C60; 92C50; 45K05; 65R20; 65C05.

1 Introduction

Infection dynamics are typically shaped by heterogeneous and interconnected factors, including exogenous individual-based aspects such as the individual response to threat, the collective compliance to non-pharmaceutical interventions, and variable endogenous characteristics of a pathogen [6, 7, 8, 12, 21, 22, 26, 37, 44]. Indeed, the infectivity of individuals may vary strongly from the onset of the infection to its end. In a purely data-oriented approach, the day of first contact with an infected individual is typically unknown or unreported and the infectivity levels may vary in time. Therefore, the trend of an epidemic naturally involves a plurality of timescales: one that is linked to the evolution of infected cases among the population and a second related with the change of infectivity related to in-host viral dynamics [25, 33, 4].

In this direction, viral load (VL) dynamics represent a significant aspect for the disease progression and transmission. The evolution of VL represents a quantitative marker for assessing viral kinetics [20, 36, 40]. Nevertheless, the reconstruction of these trajectories is a challenging problem since data are incomplete and biased toward the peak viral load, that is typically related to the onset of symptoms in infected patients. We mention in this direction the empirical studies, see e.g. [23]. Hence, VL kinetics are strongly affected by uncertainties stemming from the data assimilation processes. Existing methods are not well-suited to consider possible mutations of the virus, resulting in a modification of the infectiousness curve associated with each individual illness. Therefore, a purely data-oriented approach should be complemented by a multiscale model-oriented approach to provide robust long-term predictions of epidemiological dynamics [28]. In the last few years, several mathematical approaches have been designed to assess the impact of infectivity dynamics at the population level [14, 32]. These works are based on kinetic-type equations and on simplified dynamics for which observable quantities are analytically computable. The developed methods are capable to link agent-based dynamics to available data, providing effective transition rates in macroscopic compartmental models, see also [3, 5, 17, 19, 43, 44].

In this study, we focus on the definition of a new pipeline capable of extracting information starting from viral load kinetics through a model-driven approach that couples the multiple timescales of infection dynamics. In detail, we exploit VL kinetics obtained in a series of SARS-CoV-2 infections that are not affected by exogenous influencing factors in the dynamics, like e.g. the implementation of mass vaccination policies. In detail, we interface data-oriented viral load dynamics with a simple compartmental model with age of infection [15, 42] to describe the evolution of infected cases from SARS-CoV-2 in Italy. A new approach to evaluate the dynamics of the viral load has been developed and compared with original data obtained during the initial waves of the SARS-CoV-2 pandemic.

2 Materials and methods

2.1 Samples and SARS-CoV-2 kinetics

Since February 2020, Fondazione IRCCS Policlinico San Matteo (Pavia, Italy) has been identified as reference regional center for the diagnosis of SARS-CoV-2 infections, through the analysis of respiratory samples. In these analyses, RNA was extracted by using the MGISP-960 automated workstation and the MGI Easy Magnetic Beads Virus DNA/RNA Extraction Kit (MGI Technologies, Shenzhen, China). SARS-CoV-2 RNA was detected using the SARS-CoV-2 variants ELITe MGB kit (ELITechGroup, Puteaux, France; cat. no. RTS170ING) targeting ORF8 and RdRp gene. Reactions were carried out on the CFX96 Touch Real-time PCR detection system (BioRad, Mississauga, ON, Canada).

Individuals that resulted positive for SARS-CoV-2 RNA may leave isolation only when two consecutive negative respiratory samples occurred. Thus, a series of follow-up samples were analyzed to evaluate the structure of SARS-CoV-2 viral loads and used to monitor the evolution of SARS-CoV-2 VL kinetics. The raw data for each patient consist of the days in which respiratory samples were collected and analyzed (where 0 always represents the day of the first test) and of the kinetics Ct detected by each test (where $40\ \textrm{Ct}$ means complete negativity to the virus and thus recovery). So, for example, a patient strongly positive ($\sim 20\ \textrm{Ct}$) to the first test, that becomes barely positive ($\sim 36\ \textrm{Ct}$) after 9 days and finally negative after 7 more days, would be characterized by the two triples $(0,9,16)$ for the days and $(20,36,40)$ for the kinetics. The values of the corresponding viral loads in terms of RNA $\textrm{copies}/\textrm{ml}$ of respiratory samples may then be deduced from the kinetics Ct using the empirical conversion formula

$$1\ \frac{\textrm{copies}}{\textrm{ml}}=16.6\times 10^{-0.2814\ \times\ 1\textrm{Ct}\ +\ 11.232},$$

where $16.6$ is the dilution factor needed to report VL as RNA $\textrm{copies}/\textrm{ml}$ of respiratory samples. Laboratory data were anonymized and retrospectively analyzed. The study was conducted in accordance with the Declaration of Helsinki and no clinical information were available, except for the age and gender.

2.2 Statistical analysis

Numerical variables are described as mean $\pm$ standard deviation, categorical variables as count and percentages. The comparison of the two periods of time considered is performed in terms of t-test for quantitative variables and Fisher exact test for quantitative variables. Considering the longitudinal structure of the data, they are presented as daily averaged values. The comparison between the two periods and the investigation of the influence of age and gender on the kinetics are performed in terms of a linear mixed-effects model. The temporal evolution of VL kinetics is obtained as a Gamma distribution using the Matlab built-in fit function, using the LAR method to ensure that the extreme values have a lesser influence on the result. The significance threshold is set as $p<0.05$.

3 Results

3.1 Data description

The study analyzed a total of 233 viral load kinetics, corresponding to approximately 700 respiratory samples, collected from patients across two distinct periods. The first period goes from November to December 2020, before the introduction of variants of concern (VOCs), and included 71 patients. The second period, from January to May 2021, during the circulation of VOCs, involved 162 patients.

The subjects considered had a mean age of $42.3\pm 14.7$ years, with $55\%$ of females (120) and $45\%$ of males (113). The demographic distribution, in terms of age and gender, does not statistically differ between the two periods. Using mixed-effects models, we found that age and gender did not significantly influence viral load kinetics. Time to negativization has a significant effect, showing a decreasing pattern of viral load kinetics ($\textrm{estimate}=-482193.7$, $\textrm{p-value}<0.001$).

Additionally, a significant difference in viral load kinetics was observed between the two periods, with lower values recorded in 2021 compared to 2020 ($\textrm{estimate}=-2122031.9$, $\textrm{p-value}=0.015$).

3.2 Data-oriented viral load kinetics

To understand what the proper target shape for the theoretical function $\beta$ could be (and its evolution over the two observation periods), we initially infer it from the best fit of the individuals’ average viral load, from each one of the considered temporal windows. Based on existing works in the medical literature about the shape of the SARS-CoV-2 infectiousness [11, 37], we expect $\beta$ to be characterized by a rapid growth in the few days after the first contact with an infective individual, until the infection reaches its maximum, followed by a slower decay, up to the complete negativization of the subject after several days. The issue with the raw data is that we cannot be sure that the first molecular test of each patient corresponds to their infections’s peak. Therefore, we loose any deterministic information on both the day in which the infection reaches its maximum and the value of this maximum. To rule out part of this randomness, we reorder the infection cycle of each subject based on their first negative test, rather than their first positive one. Indeed, while it remains unclear when the actual negativization happened, we know for certain that the value of the viral load on that day must be 0. In Figure 1 we present the raw data for the viral load kinetics of each individual patient ordered by the first negativization, over the two distinct periods of time during which their infectiousness was monitored.

Figure 1: Evolution of individuals’ viral load kinetics over time, with an ordering based on the first negative test taken by the infective patients. We distinguish between the two observation periods of November–December 2020, accounting for 71 patients (left), and January–May 2021, accounting for 162 patients (right).

Starting from this reordering, we consider the daily average of the viral load with respect to the number of patients that got tested on that day. We then infer the shape of $\beta$ as the best fit of these averages, weighted by the daily patients’ count. This means that the more patients got tested on one specific day, the more importance the corresponding average viral load value will have in the fitting. The outcomes of this analysis are shown in Figure 2. On average, the infection persists for 25 days at most and its peak is reached 3 to 4 days after the first contact with an infected individual (4 days for the viral loads considered from the period November–December 2020, and 3 days for those belonging to the window January–May 2021). The red curves represent the best weighted fits for the daily averaged viral loads and have the shape of a Gamma distribution

$$\hat{\beta}(\tau)=\alpha_{0}(\tau+\alpha_{1})^{\alpha_{2}}\exp\left\{-\alpha_{3}(\tau+\alpha_{1})\right\},$$ (1)

where $\alpha_{0}=10^{7}$ is a scale factor used to capture the magnitude of individuals’ viral load and the coefficients $\alpha_{1}\in\mathbb{R}$, $\alpha_{2}$, $\alpha_{3}\geq 0$ are reported in Table 1.

Figure 2: Best fit of the patients’ viral loads depending on the time of the infection, for the two monitoring periods of November–December 2020 (left) and January–May 2021 (right). The black crosses denote all the individuals’ infection cycles ordered based on their first negativization. The blue dots represent the averages of the daily viral loads. The red curves are the best fit for these averages, weighted by the number of patients counted on each day. The shaded region gives the $80\%$ and $95\%$ prediction bounds for the whole dataset.

November–December 2020’s sample
	Estimated value	$95\%$ CI
$\alpha_{1}$	0.4189	[-0.4865, 1.3245]
$\alpha_{2}$	2.2015	[1.6124, 2.7906]
$\alpha_{3}$	0.5701	[0.3928, 0.7475]

January–May 2021’s sample
	Estimated value	$95\%$ CI
$\alpha_{1}$	0.6262	[0.1588, 1.0936]
$\alpha_{2}$	1.5687	[1.0462, 2.0912]
$\alpha_{3}$	0.4964	[0.3348, 0.6580]

Table 1: Optimised coefficients for viral load kinetics, for the two observation periods of November–December 2020 (left) and January–May 2021 (right).

3.3 Epidemic model with uncertain quantities

Compartmental models are classically defined to describe mathematically the spread of epidemics in a population [16]. The main idea of such models is to divide the population into different subgroups, each identified by a specific epidemiologically relevant status. Suitable transition rates are introduced to describe the switches between compartments. Generalizations of the mentioned compartmental models include additional factors characterizing the transmission dynamics, such as external influence, age structure, variable contact dynamics and mobility of agents, see e.g. [1, 2, 9, 18]. In the following, we will consider a SIR-type model with age of infection where the infectivity keeps track of the disease’s variability inside each individual and evolves in parallel with the main timescale of the epidemic. This allows to introduce a local incidence rate function of the age of infection, modelling the average viral load of the infected individuals and whose shape can be inferred by the information coming from real data on molecular tests. To take into account the structural uncertainties stemming from data-assimilation processes of VL kinetics, we extend the model approach to consider the presence of uncertain quantities as a structural feature of the epidemic dynamics. To date, few results are available regarding the development of uncertainty quantification (UQ) methods in epidemic systems, see [1, 10, 13, 38]. The development of UQ methods are based on an increased dimensionality of the problem which is made dependent on a random variable $\mathbf{z}\in\mathbb{R}^{d_{z}}$ of which the distribution $p(\mathbf{z})$ is known [24, 34, 41]. The extrapolation of statistics is classically obtained by looking at quantities of interest, with respect to $\mathbf{z}$, representing the expected solution of the problem.

At a given time $t\geq 0$, we subdivide the population into susceptible $S$, recovered $R$ and infected with variable infectiousness $i$. We then consider the new triple $(S(t,\mathbf{z}),i(t,\tau,\mathbf{z}),R(t,\mathbf{z}))$ of unknowns, solution to

$$\left\{\begin{array}[]{l}\displaystyle\frac{\mathrm{d}S(t,\mathbf{z})}{\mathrm{d}t}=-\bar{\beta}S(t,\mathbf{z})\int_{\mathbb{R}_{+}}\beta(\tau,\mathbf{z})i(t,\tau,\mathbf{z})\mathrm{d}\tau,\\[14.22636pt] \displaystyle\partial_{t}i(t,\tau,\mathbf{z})+\partial_{\tau}i(t,\tau,\mathbf{z})=-\bar{\gamma}\gamma(\tau,\mathbf{z})i(t,\tau,\mathbf{z}),\\[14.22636pt] \displaystyle\frac{\mathrm{d}R(t,\mathbf{z})}{\mathrm{d}t}=\bar{\gamma}\int_{\mathbb{R}_{+}}\gamma(\tau,\mathbf{z})i(t,\tau,\mathbf{z})\mathrm{d}\tau,\end{array}\right.$$ (2)

where now $\beta(\tau,\mathbf{z})$ and $\gamma(\tau,\mathbf{z})$ represent the uncertain infectiousness and recovery of the infected individuals, depending on the age of infection $\tau\geq 0$, while the constants $\bar{\beta}>0$ and $\bar{\gamma}>0$ enclose the average contact rate and duration of the illness. We observe that the observable number of infected also varies based on the new timescale $\tau\geq 0$, meaning that the quantity $i(t,\tau,\mathbf{z})\geq 0$ depends on both the timescale, and it is related to the total number of infected $I$ through the relation

$$I(t,\mathbf{z})=\int_{\mathbb{R}_{+}}i(t,\tau,\mathbf{z})\mathrm{d}\tau.$$

This model encapsulates the timescale of the epidemic with the one characterizing the evolution of the infection, in terms of infectiousness. For each $\mathbf{z}$, the system is finally completed with suitable initial and boundary conditions for the variables $S(t,\mathbf{z})$, $i(t,\tau,\mathbf{z})$ and $R(t,\mathbf{z})$. A natural choice is

$$S(0)=S_{0},\quad i(0,\tau,\mathbf{z})=i_{0}(\tau,\mathbf{z}),\quad R(0)=R_{0},$$ (3)

$$i(t,0,\mathbf{z})=\bar{\beta}S(t,\mathbf{z})\int_{\mathbb{R}_{+}}\beta(\tau,\mathbf{z})i(t,\tau,\mathbf{z})\mathrm{d}\tau,$$ (4)

where $S_{0}$, $i_{0}(\tau,\mathbf{z})$, $R_{0}\in\mathbb{R}_{+}$ are such that

$$S_{0}+\int_{\mathbb{R}^{d_{z}}}\int_{\mathbb{R}_{+}}i_{0}(\tau,\mathbf{z})\mathrm{d}\tau\mathrm{d}p(\mathbf{z})+R_{0}=1,$$

where $p(\mathbf{z})$ is the distribution of the random variable $\mathbf{z}$.

3.4 Inference of viral load dynamics based on uncertain data

In this section we define a pipeline to effectively estimate actual viral load kinetics $\beta(\tau,\mathbf{z})\geq 0$, by using available epidemiological data and introducing uncertainty in the day of first contact, expressed by the coefficient $\alpha_{1}$ in (1). To this end, we consider an uncertain parameter $\alpha_{1}(\mathbf{z})\geq 0$ to make the local incidence rate random. Moreover, we take into account data from the COVID-19 pandemic, both at local and national level, by considering the evolution of infected in the province of Pavia and then extending our analysis to the whole territory of Italy, where the central government implemented several non-pharmaceutical interventions (NPIs) during the pandemic. In Table 2 we report the timeline of the main NPIs taken before the mass vaccination campaign started in mid 2021. It is however worth mentioning that local additional measures, like isolation of small portion of a territory, had often been implemented as well.

Before advancing to the estimation of VL kinetics, it is crucial to address some key points. We stress the inherent complexity in calibrating epidemiological models, which is particularly demanding from a numerical perspective. Indeed, the process of aligning these models with real-world scenarios requires a careful approach to parameter estimation. Furthermore, available data on viral load kinetics tend to hide the initial evolution of this dynamics, since agents are typically tested after the manifestation of the symptoms. This bias is primarily attributed to limited testing capacity, resulting in a dataset that often reflects a conservative estimate. Consequently, any attempts to estimate the viral loads should be approached with a mindful consideration of these inherent limitations, coupled with a commitment to refining the model’s accuracy through continuous reassessment and adaptation.

Implemented NPI	Date	Average recovery rate
First lockdown	March 9, 2020
Curfew measures	October 22, 2020	14 days
Strengthening of measures	March 10, 2021

Table 2: Timeline of the main NPIs employed in Italy before the mass vaccination campaign started in mid 2021.

We shall focus on three different periods of time. We start by analyzing the first COVID-19 epidemic wave in Italy, which began in February–March 2020. We then proceed with the study of the second wave of October 2020, and conclude with the third major wave of February–March 2021. In particular, we point out that our epidemiological data for the province of Pavia were only acquired up to January 18 2021, and we will therefore narrow down our analysis on the third wave to the sole case of Italy, for which we shall instead take advantage of the open access data collected by the Johns Hopkins University (github.com/CSSEGISandData/COVID-19), reporting in details the daily number of infected, recovered and deceased individuals in Italy up to August 4 2021. Let $\mathcal{C}_{z}=\{\mathbf{z}_{1},\dots,\mathbf{z}_{M}\}$ be a sample of the random variable $\mathbf{z}$ with probability density function $p(\mathbf{z})$. We have adopted a multi-level approach for estimating the parameters $(\alpha_{i})_{i=1,2,3}$ that characterize $\beta(\tau,\mathbf{z}_{j})\geq 0$, which is articulated in two steps.

Step 1 – Determination of the epidemiological parameters. For starting, we initially assess the (deterministic) value of $\alpha_{1}$ dictated by the epidemic during the two weeks preceding any of the three main NPIs presented in Table 2, and we estimate at the same time the relevant epidemiological parameter $\bar{\beta}\geq 0$ which defines the baseline transmission dynamics, assuming no restrictions on social contacts were in place. Notice that we have to discard the variable $\mathbf{z}$ here, because we first need to ensure that the day of first contact $\alpha_{1}$ is compatible with the epidemiological data, in order to quantify the shift between the infection and the epidemic timescales. For this, we look for the optimal $\alpha_{1}$ and $\bar{\beta}$ minimizing the relative weighted $L^{2}$ norm between the reported number of infected and recovered individuals (in Pavia or in Italy, depending on the local or national setting under consideration) which we denote $\hat{I}(t)$ and $\hat{R}(t)$, and the theoretical evolution of these quantities, $I(t)=\int_{\mathbb{R}_{+}}i(\tau)\mathrm{d}\tau$ and $R(t)$, given by the model (2) in absence of uncertainty. More precisely, we solve the following minimization problem

$$\min_{\alpha_{1},\bar{\beta}}\ (1-\eta)\frac{\|I(t)-\hat{I}(t)\|_{L^{2}([t_{0},t_{f}])}}{\|\hat{I}(t)\|_{L^{2}([t_{0},t_{f}])}}+\eta\frac{\|I(t)+R(t)-\hat{I}(t)-\hat{R}(t)\|_{L^{2}([t_{0},t_{f}])}}{\|\hat{I}(t)+\hat{R}(t)\|_{L^{2}([t_{0},t_{f}])}},$$ (5)

over the time horizon $[t_{0},t_{f}]$ representing any of the three periods that can be identified with the initial spreading phase of the three epidemic waves under study: February 24–March 9 2020, October 7–October 22 2020, and February 24–March 10 2020. In particular, we choose here the infectiousness $\beta(\tau)$ to be the Gamma functions (1), obtained in the previous section from the observed viral loads. We use the fitting from November–December 2020 as initial guess for the optimization (5) during the first two epidemic waves, while we opt for the January–May 2021’s fitting to study the third wave. Notice that at this stage the parameters $\alpha_{2}$ and $\alpha_{3}$ are fixed from these two fits, since we only need to ensure a proper initial calibration of the time shift $\alpha_{1}$, to correctly align our model with the epidemiological data. Finally, we always consider a constant recovery rate $\bar{\gamma}=1/14$ to lower the complexity of the minimization. This value is consistent with the existing literature on the COVID-19 infection [12, 30, 44], for which the time to viral clearance has been reported to span approximately between 10 and 21 days on average over the different variants. In Table 3, we collect all the parameters obtained from this initial analysis, for each of the three waves.

Province of Pavia
	First wave	Second wave
$\alpha_{1}$	1.8110	-0.6706
$\alpha_{2}$	2.2015	2.2015
$\alpha_{3}$	0.5701	0.5701
$\bar{\beta}$	0.1486	0.1603
$\bar{\gamma}$	1/14	1/14

Italy
	First wave	Second wave	Third wave
$\alpha_{1}$	6.6849	-0.6032	-0.8251
$\alpha_{2}$	2.2015	2.2015	1.5687
$\alpha_{3}$	0.5701	0.5701	0.4964
$\bar{\beta}$	0.3992	0.1225	0.1323
$\bar{\gamma}$	1/14	1/14	1/14

Table 3: Epidemiological parameters and coefficients of the infectiousness fuction, for the three periods of time under consideration, both at a local (province of Pavia) and a national level (Italy).

Step 2 – Inclusion of the uncertainty. We now go further in the analysis by including the presence of uncertainty, by looking at the full model (2). The aim is then to perform an uncertainty quantification on the shape of $\beta(\tau,\mathbf{z}_{j})$, with respect to the randomness $(\mathbf{z}_{j})_{j=1}^{M}$ that may appear when estimating $\alpha_{1}(\mathbf{z}_{j})$, the day of first contact. For this, we perturb the optimal value $\alpha_{1}$ computed in the previous step by a centered uniform random distribution $\mathcal{U}([-2,2])$, i.e. $\alpha_{1}(\mathbf{z}_{j})=\alpha_{1}+\mathbf{z}_{j}$ with $\mathbf{z}_{j}\sim\mathcal{U}([-2,2])$, and we assess the impact of this addition on the parameters $\alpha_{2}(\mathbf{z}_{j})$ and $\alpha_{3}(\mathbf{z}_{j})$ that define the shape of the $\beta(\tau,\mathbf{z}_{j})$. The (randomized) optimal shape is thus determined by fitting model (2) to reproduce the data on infected and recovered cases, through a minimization problem similar to (5), but with parameters $\alpha_{1}(\mathbf{z}_{j})$ and $\alpha_{2}(\mathbf{z}_{j})$ instead of $\alpha_{1}$ and $\bar{\beta}$, while of course the latter and $\bar{\gamma}$ are those obtained from the previous step. This procedure is repeated for the whole sample $\mathcal{C}_{z}$, and we can finally infer the optimal shape of the infectiousness function with UQ, to be the expected value $\mathbb{E}_{z}(\beta(\tau,\mathbf{z}))$ with respect to $\mathbf{z}$ and based on the probability density $p(\mathbf{z})=\mathcal{U}([-2,2])(\mathbf{z})$.

Figure 3: Data- and model-driven inference of the optimal shape of the viral load in the province of Pavia (Italy), over two different epidemic waves: the one occurred at the beginning of the pandemic in 2020, with a focus on the two weeks (February 24–March 9) before the first lockdown was established (top), and the other occurred in the fall of the same year, with a focus on the two weeks (October 7–22) preceding the imposition of a national curfew (bottom). The figures on the left provide the optimal viral load shapes (red curves) determined for these two periods, starting from the data collected by Fondazione IRCCS Policlinico San Matteo. The shaded areas show how the infectiousness function is altered when a uniform noise $\mathcal{U}([-2,2])$ modifies the day of first contact. The figures on the right present the evolution of infected from SARS-CoV-2 in the province of Pavia over the considered windows of time. We compare computed (in red) and reported (in black) number of cases.

The outcomes of our two-step procedure are reported in Figure 3 for the case of Pavia and in Figure 4 for the case of Italy. In the left figures we plot in red the curves providing the mean best shape of the infectiousness function over the three epidemic waves that we have considered. The shaded regions show the $80\%$ and $95\%$ percentiles of the different shapes, based on the sample $\mathcal{C}_{z}$. We then perform numerical simulations to assess the viability of the model (2) to predict the evolution of infected individuals during these subsequent phases of the pandemic. We take the infectiousness function to be the computed mean best shapes over the three periods and we choose $\gamma(\tau,\mathbf{z})$ constant and equal to 1 for simplicity. In particular, the duration of the infection is solely determined by $\bar{\gamma}$ and the contact rate is given by the optimal ones $\bar{\beta}$ determined from Step 1 and reported in Table 3. We provide in Figures 3 and 4, on the right, the comparisons between the reported number of infected in Pavia and in Italy over the three epidemic waves, and the corresponding predictions of the model (2), with the aforementioned choices. We notice that with the progression of the pandemic, and specifically over the third wave in March 2021, the model loses some accuracy since it is less appropriate to describe slower growths of the epidemic, stemming from the numerous containment measures introduced by the Italian government in the course of 2021 to reduce contacts among the population. This suggests in particular that in order to capture such intricate effects and minimize the error made by the model, one should consider a variable contact rate $\bar{\beta}(t)$ that continuously evolves over the epidemic time $t$.

Figure 4: Investigations on the SARS-CoV-2 pandemic in Italy, during three distinct epidemic waves. Top row: analysis of the epidemic wave that started in February–March 2020. Optimal shape $\mathbb{E}_{z}(\beta(\tau,\mathbf{z}))$ of the infectiousness function (left) and model predictions (right) for the cases of infection registered during the period February 24–March 9. Middle row: analysis of the new epidemic outburst that took place in October 2020. Optimal shape of the viral load (left) and model predictions (right) for the cases of infection registered during the period October 7–22. Bottom row: analysis of the epidemic wave occurred from March 2021. Optimal shape of the viral load (left) and model predictions (right) for the cases of infection registered during the period February 24–March 10.

Extending the model predictions over the whole timeframe. We now perform a calibration of our model to cover the whole timespan of observation, connecting the epidemic waves between them and providing full predictions for the evolution of infected. We initially observe that as the epidemiological parameter $\bar{\beta}$ models the contact rate between individuals, its value is highly influenced by the restrictions imposed by a government. This fact needs to be taken into account when extending our model’s predictions from one epidemic wave to the other, since various containment measures (we mention in particular the lockdown imposed between March and May 2020) have been put in place by the Italian government over time to slow down the epidemic spreading. Since we do not need to evaluate the shape of the viral load in these connecting timeframes (the previous two-step analysis is relevant only when the virus is able to spread free among the population, allowing for a proper calibration of the epidemiological parameters), in order to connect the observations we can rely on a simpler model that allows to speed up the computational time. Specifically, we make use of system (2) with the choices $\beta(\tau,\mathbf{z}),\ \gamma(\tau,\mathbf{z})\equiv 1$. From the above considerations, we account for the variable contact restrictions by assuming that the contact rate varies over the epidemic time as $\bar{\beta}=\bar{\beta}(t)$. Then, we determine the value $\bar{\beta}(t)$ by solving an optimization problem for a sequence of time steps $t_{i}$ (days) over a moving time window of three days (namely, averaging the fitting over three days), varying in any timespan $t\in[t_{0}+1,t_{f}]$ where $t_{0}$ and $t_{f}$ correspond respectively to the end of an epidemic wave and the beginning of the next one. This strategy is needed in order to correctly assess the evolution of the contact rate depending on the containment measures that are in place. More precisely, we solve the following constrained least-square problem

$$\min_{\bar{\beta}(t_{i})}\ (1-\eta)\frac{\|I(t)-\hat{I}(t)\|_{L^{2}([t_{i}-1,t_{i}+2])}}{\|\hat{I}(t)\|_{L^{2}([t_{i}-1,t_{i}+2])}}+\eta\frac{\|I(t)+R(t)-\hat{I}(t)-\hat{R}(t)\|_{L^{2}([t_{i}-1,t_{i}+2])}}{\|\hat{I}(t)+\hat{R}(t)\|_{L^{2}([t_{i}-1,t_{i}+2])}},$$ (6)

and the recovery rate $\bar{\gamma}$ is always kept fixed to the value $1/14$. We complete our study by pushing the model predictions up to the last observation times of our experimental data on the epidemic, that are January 18 2021 for the case of the province of Pavia and August 4 2021 for the case of Italy. Using the above procedure, we connect the observations by calibrating $\bar{\beta}(t)$ with the minimization strategy (6). Using all these optimized parameters, we can finally run our model over the whole timeframe of observation. We depict the resulting predictions for the evolution of infected individuals in the province of Pavia, in Figure 5 and Figure 6 respectively, on the right. The confirmed cases are plotted in black, while the model approximations are shown in red.

Figure 5: Left: comparison between viral loads from the first (in orange, March 2020) and the second (in blue, October 2020) epidemic waves, in the province of Pavia. In order for the SIR model (2) to best fit the data on confirmed cases, the peak of the infection should be reached about 6 to 8 days after the first contact with an infected subject. Right: evolution of infected from SARS-CoV-2 in Pavia from February 24 2020 to January 18 2021. Comparison between computed (in red) and reported (in black) number of cases.

At last, in Figures 5 and 6 on the left, we provide a direct comparison between the optimal viral loads from different epidemic waves, obtained via the previous two-step analysis, together with their $80\%$ prediction bounds. These are aligned by their respective peak of infectiousness (which corresponds to day 0 on the graph). One can observe that the infection’s peak progressively decreases, while the tails of the viral load tend to grow. This suggests that the virus becomes weaker but more effective over a longer period of time, and individuals stay infective for more days on average. Moreover, the decrease in curve’s growth in the beginning of the infection could suggest a slower progression of the symptoms and explain the greater spreading of the epidemic with each new wave.

Figure 6: Left: comparison between viral loads from three distinct epidemic waves in Italy, happened respectively in March 2020 (orange), in October 2020 (blue) and in March 2021 (green). In order for the SIR model (2) to best fit the data on confirmed cases, the peak of the infection should be reached about 5 to 8 days after the first contact with an infected subject. Right: evolution of infected from SARS-CoV-2 in Italy from February 24 2020 to August 4 2021. Comparison between computed (in red) and reported (in black) number of cases.

4 Discussion

Understanding the structure and evolution of VL kinetics is crucial in the context of public health, when setting up quarantine and isolation measures of positive subjects. In particular, one of the critical issues emerged since the onset of the recent SARS-CoV-2 pandemic has been the need of criteria to assess the infectivity of patients and decide their releasing from prolonged quarantine, safely admitting them back to work duties and social activities. It was demonstrated that high Ct values (low VL), greater than 30-35 Ct, were associated with a reduced infectivity of the patients [29, 39, 35]. Nevertheless, existing studies indicate that the peak viral load is linked closely with the onset of symptoms making the available datasets partial and of difficult interpretation. In this work, we have established a pipeline to assess the evolution of in-host VL from respiratory samples of patients. To this end, we combined the use of real data and model-based methods. Furthermore, the provided method is capable of tracing the whole VL kinetics. In particular, our model helps to provide a data-oriented understanding of the non-infectious status of people that were under quarantine with high Ct.

Samples included in this study were collected over two periods of time (November–December 2020 and January–May 2021), that include the initial circulation of VOCs and before the extensive introduction of vaccination. We have first determined the shape of the average VL kinetics hinted by these data (recovering Gamma-type distributions) and we have successively introduced a suitable compartmental model of SIR-type, using a local incidence rate function that varies with the age of infection, to account for the evolution of the epidemic alongside that of the virus. The system of equations belongs to a large class of epidemiological models detailing the transmission dynamics through the inclusion of additional factors such as external influence, age structure, variable contact dynamics and mobility of agents [1, 2, 9, 18, 27]. In this work, we exploit the information coming from real data on molecular tests and on the reported number of infected individuals to deduce the optimal shape of the infectivity function, based on a combined data-driven and model-driven approach. To take into account the extreme variability of VL dynamics, depending on disease severity and patients’ features as well [31, 44], we also make use of UQ methods [24, 34, 41] adapted for the study of epidemic systems [1, 10, 13, 38], to take into account the uncertainties inherent to the data-assimilation processes of VL kinetics. The model allows us to infer a plausible evolution of the SARS-CoV-2 infectiousness over three distinct periods of the pandemic (February–March 2020, October 2020 and February–March 2021), in agreement with existing works from the medical literature [11, 37]. In particular, we observed that the peak of the infection function decreases over time while its tails tend to grow, suggesting that the virus becomes progressively weaker but more effective, since individuals remain infective for more days on average. This change in shape strongly manifests during the third observation period (February–March 2021) when the Alpha variant was actively circulating in Italy, providing an explanation for the greater spreading of the contagion at the time.

In conclusion, our work aimed to implement strategies of preparedness for a new pandemic, since the model and methods developed here could be easily adapted to investigate real data coming from epidemics that are associated with other respiratory viruses (e.g influenza and respiratory syncytial virus), in order to assess the evolution of the associated viral infectiousness that may lead to more informed decisions on preventive measures to reduce the spreading of the disease.