A mathematical model for the within-host (re)infection dynamics of SARS-CoV-2

Interactions between SARS-CoV-2 and the immune system during infection are complex. However, understanding the within-host SARS-CoV-2 dynamics is of enormous importance, especially when it comes to assessing treatment options. Mathematical models have been developed to describe the within-host SARS-CoV-2 dynamics and to dissect the mechanisms underlying COVID-19 pathogenesis. Current mathematical models focus on the acute infection phase, thereby ignoring important post-acute infection effects. We present a mathematical model, which not only describes the SARS-CoV-2 infection dynamics during the acute infection phase, but also reflects the recovery of the number of susceptible epithelial cells to an initial pre-infection homeostatic level, shows clearance of the infection within the individual, immune waning, and the formation of long-term immune response levels after infection. Moreover, the model accommodates reinfection events assuming a new virus variant with either increased infectivity or immune escape. Together, the model provides an improved reflection of the SARS-CoV-2 infection dynamics within humans, particularly important when using mathematical models to develop or optimize treatment options.


Introduction
The interactions between severe acute respiratory syndrome coronavirus type 2 (SARS-CoV-2) and the immune system during infection are complex.Viral load measurements of the upper respiratory tract provide a quantitative method to capture the within-host SARS-CoV-2 infection dynamics [33].Mathematical models have been developed to dissect and quantify the underlying pathogenic mechanisms, incorporating features such as viral loads, time of clinical symptom onset or infectiousness [7].
Viral within-host dynamics are commonly described by the target-cell limited (TCL) model [2].The TCL model describes the dynamics of susceptible (target) cells, infected cells and free virus populations during the course of an infection.In the past few years, TCL-based models have been employed to also describe the within-host infection dynamics of SARS-CoV-2.Early efforts focused on quantitatively describing the measured within-host viral load dynamics in individuals [3,32] and on linking viral loads to infectiousness [15,21].These models differentiated between infectious and non-infectious virus, such as residual viral genome fragments, or by adding simple, often nonmechanistic descriptions of the human immune response [3,15,21,32].By adding an interferon mediated refractory cell population or reduction in the infection and viral production rates, later models integrated more mechanistic formulations of the human immune response [14,20].These models provided further insights into SARS-CoV-2 infection dynamics, revealing heterogeneity in the infectiousness levels between individuals [14] and suggesting differences in viral dynamics between variants [20,22].TCL-based models of SARS-CoV-2 infection dynamics have lately also been used to describe the viral load dynamics during antiviral treatment and, consequently, to obtain insights into phenomena such as viral rebound and resistance formation [24,25].To date, models describing the within-host SARS-CoV-2 infection dynamics typically focus on the short-term acute phase of an infection.The assumptions that these models make are often at odds with simple but highly relevant and important clinical realities, such as the re-establishment of susceptible cells to pre-infection levels, full and permanent clearance of the infection from the individual, the corresponding decline of the immune response, and the associated established long-term post-acute infection immunity levels after exposure and infection.Here, we present a within-host SARS-CoV-2 model, which offers a more realistic account of the underlying infection dynamics.While our model captures the short-term within-host infection dynamics just as well as existing models and demonstrates a good generalizability, it considerably improves earlier work by further capturing expected post-acute infection dynamics.These include the re-establishment of the number of susceptible epithelial cells to the initial pre-infection level, a full and permanent clearance of the infection from the individual, the formation of long-term postacute infection immune response levels after infection, and the waning of immunity.Finally, we demonstrate that the model can accommodate reinfection events, taking into account differences in infectivity and levels of immune escape of the reinfecting virus variant.

Within-host SARS-CoV-2 model
Our model comprises four quantities corresponding to the populations of susceptible epithelial cells (S) and infected epithelial cells (I), free virus particles (V ), and the immune response (B, normalized to a range between 0 and 1) (Figure 1A).In the case of reinfection with a different variant, the immune response is re-activated.For a different variant with immune escape (yellow), the effectiveness of the initial immune response against the different variant suffers a drop due to the variant's escape properties.For a different variant featuring higher infectivity (blue), the immune response rises to higher levels in order to successfully control the infection (blue).The dotted line represents the immune response without reinfection.
The full system of non-linear ordinary differential equations (ODEs) for the susceptible (S) and infected (I) cells, free virus (V ), and the immune response (B) is given by the following equations: where r(t) = 1 − B(t) is the immune response effect function, and S(0) = S 0 , I(0) = I 0 , V (0) = V 0 , and B(0) = 0 are the given initial values.In a disease-free system, susceptible cells are in homeostasis, where cells are constantly produced at rate p S and die at a natural death rate d S , where p S = S 0 × d S .Once a free virus is introduced, the system is perturbed as follows: When in contact with free virus, susceptible cells get infected according to an overall infectivity rate β 0 r(t), which is the product of the infectivity rate constant β 0 and immune response effect function r(t).
The immune response effect function indicates the remaining proportion of successful infections of susceptible cells despite the presence of the immune response.The term of the first equation in ( 1) is to be interpreted as (i) the removal of susceptible cells due to infection assuming no activity of the immune response (B(t) = 0) and (ii) the immediate re-introduction of susceptible cells involved in unsuccessful infections inhibited by the anti-SARS-CoV-2 activity of the immune response.The successfully infected cells enter the infected cell population and are eliminated from that population with virus-induced death rate d I .Due to the non-lytic property of SARS-CoV-2, infected cells constantly produce and release free virus at rate p V , which is cleared at rate d V [8].Free virus leading to successfully infected cells are lost from this population according to the overall infectivity rate β 0 r(t).A SARS-CoV-2 infection activates the immune response (Figure 1B).First, the innate immune response detects the viral infection and limits its spreading within the host.The innate immune response also triggers the adaptive immune mechanisms, which ultimately clear the virus and form a lasting, long-term SARS-CoV-2-specific immune memory [4,27,30].This ensures an enhanced response to future infections of the same virus.The term immune response B in our model represents the combined effects of innate and adaptive immunity on the ability to control the virus.For a primary infection, the individual was assumed to not have been exposed to SARS-CoV-2 before, and hence, the immune response is initially described by a naive pre-infection state with B(0) = 0. Upon infection, B increases proportional to the viral load at rate p B and is limited by the upper bound of 1, such that B(t) ∈ [0, 1].The activation of the immune response is assumed to lead to a reduced overall infectivity and hence, to fewer successful infections of susceptible cells [15,20].Long-term, we assumed that the immune response wanes to a low level of immunity after infection.To counteract a repeated spreading of the same virus within an individual, we defined a minimal immune threshold B thres against which the immune response B converges to in the longterm with rate d B , while assuring a permanent clearance of the infection.The immune threshold B thres is given by: How B thres is determined is described in detail in the section Methods.To keep the immune response inactive during the initial virus-free state and allow for its activation upon contact with SARS-CoV-2 only, we multiplied d B (B(t) − B thres ) by B(t) in the last equation of (1).In summary, the overall infectivity rate is  2A).The main differences between model fits are in the predicted initial infection dynamics prior to peak viral loads.From the estimated individual-specific viral production rate constants p V and the fixed death rate constant of an infected cell d I , we calculated an estimated median number of 75 viral particles (interquartile range [66, 98]) being produced during the life span of an infected cell.The average produced number of viral particles per infected cell and individual is given by p V d I .Furthermore, we qualitatively compared the fits of our model to those of Ke et al. for all 56 individuals, using the mean squared errors (MSEs) of the respective model fits (Figure 2B).The MSEs are closely scattered around the diagonal line, suggesting a good general agreement between the two models across individuals.For 33 of the 56 individuals, our model fits led to lower MSEs than those of Ke et al. (dots in orange area of Figure 2B).Furthermore, we considered the Bland-Altman method for statistically assessing the agreement between the model fits [13].Considering the log 10 MSEs, we found a negligible bias (log 10 MSE -log 10 MSE Ke ) and identified that in 95% of the individuals the MSEs of Ke et al. are between 0.46 and 1.85 times the MSEs of our model (Figure 2C).Overall, the two models concur well in capturing the SARS-CoV-2 infection dynamics of nasal viral load samples in humans.

Numerical results
To test the generalizability of the model, we performed a simulation study.Similar to the model fitting in the previous section, we assumed the rate constants p V , p B , and d B to be individual-specific and to capture the heterogeneity underlying the SARS-CoV-2 infection dynamics.All other rate constants were assumed to be fixed and to be the same across individuals.We drew 50 rate constant triplets out of the multivariate distribution we received from the estimated individual-specific rate constants in the previous section.Using these sets of rate constants and our within-host model, we simulated 50 SARS-CoV-2 infection dynamics (Figure 3A).From all 42 fits and 50 simulations, we computed five key features describing the underlying infection dynamics, namely (i) the log 10 peak viral load (VL), (ii) B thres , (iii) the number of days from infection to peak VL, (iv) the number of days from peak VL to undetectable VL (defined as a CN value ¿ 42), and (v) the number of days from infection to B ≥ B thres .We found the distributions of each of these five key features to be statistically comparable, that is we cannot reject the null hypothesis that the medians of both groups are equal, between the fitted and simulated infection dynamics with p-values 0.37, 0.51, 0.42, 0.52, and 0.45, respectively (Figure 3B).The median values for each of the key features for our fitted and simulated infection dynamics are shown in Table 1  Our model and simulation output agree well with those reported by others.Together, the model reliably generates simulated individual-specific SARS-CoV-2 infection dynamics.

Post-acute infection and reinfection dynamics
To verify the model performance for biological realism of its output, we further considered our model fits of the nasal viral load data (Ke et al. [14]) for a period of 90 days post infection (Figure 4).For all but one individual, the susceptible cells S re-establish their initial pre-infection levels after infection, which is denoted by S = 1 at 90 days since infection (Figure 4, left panel).Moreover, the infection is permanently and fully cleared within all individuals, given by I = 0 and V = 0 after 90 days since infection (Figure 4, middle panels).Finally, the model fits of all individuals account well for immune waning and the establishment of long-term post-acute infection immune response levels at B thres after infection (Figure 4, right panel, individual-specific B thres -values not shown).
Through its intrinsic properties, our within-host SARS-CoV-2 model captures effectively the expected post-acute infection dynamics with respect to the different model compartments.This is an improvement with respect to earlier models, which provide anomalous predicted infection dynamics when extended to a long-term post-acute infection period.According to the USA Centers of Disease Control and Prevention, a reinfection is a positive detection of SARS-CoV-2 at least 90 days after a previous infection [11].At present, reinfections typically involve a different SARS-CoV-2 variant and hence, potentially different infection dynamics due to differences in both the reinfecting variant or the variant-specific immune response [23].In our model parameterization, we accounted for both of these factors.We explored a scenario, where the new variant features an increased infectivity, which is reflected by an increased infectivity rate constant β 0 in our model (Figure 5A, left panel).Furthermore, we explored a scenario, where the new variant features immune escape.In our model, immune escape is reflected by an initial drop in the variant specific immune response B (Figure 5B, left panel).We used the previously simulated post-acute primary infection dynamics of a randomly selected representative individual from the Ke et al. data set.Furthermore, we assumed the individual to have been in contact with a different variant with either increased infectivity or immune escape, leading in each of these scenarios to a single successfully infected cell at 90 days since primary infection.Using the adapted variantspecific model parameterizations, that is an increased infectivity rate constant of β0 = 2β 0 or a decreased initial immune response of B(0) = 0.5B(90), where β0 the variant-specific infectivity rate constant and B(t) the variant-specific immune response, we simulated the viral loads and immune responses for both variants for another 60 days after reinfection (Figure 5A and B, middle and left panels).Our corresponding simulation results show a detectable increase in the viral loads due to reinfection for both variants.As a result of pre-existing immunity at the time of reinfection and given the parameterizations of β0 and B(0), the peak viral loads remain lower during reinfection than observed during primary infection for both variants.Overall, the model accommodates well for reinfection events and differentiates between reinfection modes involving distinct variant properties.

Discussion
We developed a mathematical model that describes the within-host SARS-CoV-2 infection and reinfection dynamics.Our model efficiently captures the acute short-term viral infection dynamics of the virus, and importantly, it further successfully accommodates biological phenomena surrounding the period following acute infection.By including viral loads and immune response dynamics during reinfection events by SARS-CoV-2 variants with increased infectivity or immune escape properties, our model accounts for aspects of the within-host infection process that go beyond current approaches.We fitted the model to data from a previous clinical study that measured viral loads of infected individuals and found the model to describe the acute short-term infection dynamics equally well as the model fits of Ke et al. (Figure 2).We fixed all model parameters for which there are reliable value estimates available from the literature and only estimated the three rate constants for which literature yields a very wide range of possible values or no values at all.These three rate constants are the viral production rate constant and the two rates determining the dynamics of the immune response: the activation and the waning rate constants.Numerical results showed that the fits of our model to the individual viral load are sensitive to these three model parameters.Small differences in the initial infection dynamics prior to peak viral loads may be the result of our model assuming that peak viral load occurs at 6 days post infection, while Ke et al. estimated the duration from infection to peak viral load.As there is only little data collected before the peak viral load, it is difficult to comment on these discrepancies between models.Using our model fits, we estimated that a single infected cell produces on average 75 viral particles in total (interquartile range [66, 98]) during its life span.The number of infectious viral particles cells produce is an important property in studying intra-host viral dynamics.Our value is in the range of those reported by earlier studies: from 10 to 100 infectious particles produced by a single infected cell [9,29].Due to the time-dependent description of the immune response and the introduction of an immune threshold, B thres , which is the minimal immune response required to counter a repeated spreading of the same virus variant within the individual, our approach provides insights into additional features regarding the course of infection within an individual.In the post-acute infection phase, the model recapitulates the restoration of the susceptible cell numbers to pre-infection homeostatic levels and accounts for the complete clearance of the infection, for which the immune threshold is critical (Figure 4).After an acute infection, the immune response B converges to the immune threshold and, hence, maintains at long-term post-acute infection levels.The immune threshold in our model may reflect the persistence of the underlying immune response across sequential SARS-CoV-2 infections as observed by Kissler et al. [17].For the definition of B thres , we determined the basic reproduction number R 0 , when no immune response is present, and an effective reproduction number R that is immune response-dependent and changes over time (Figure 6).Stability analysis of the model allows for a lucid description of the dynamic behaviour of the populations involved in the infection process and their steady states.Together, through its intrinsic properties, our within-host SARS-CoV-2 model captures effectively the expected long-term post-acute infection dynamics with respect to the different model populations.This is an improvement on earlier models, which when extended to the post-acute infection period provide anomalous predicted infection dynamics.
Our model further incorporates reinfection events by SARS-CoV-2 variants with different properties affecting within-host viral dynamics (Figure 5).Introducing a new virus variant with increased infectivity or immune escape perturbs the convergence of B to the immune threshold by raising the threshold or by lowering the effective immunity, respectively.As a result, the immune response against the new virus variant at the time of secondary infection is lower than the new immune threshold (Figure 4).According to the formal mathematical definition of the immune threshold of the virus variant Bthres in our model, the virus variant spreads if B < Bthres (Figure 6).Only when B reaches a value greater than Bthres is the spread of the new virus variant contained by the new updated immune response B. For both kinds of new virus variants, this feature allows the model to account for the initiation of the reinfection.That ability of our model to accommodate reinfection events and differentiate between reinfection modes involving distinct variant properties allows for its applications in clinical and public health scenarios, providing a more realistic description of key outcomes.This study has several limitations.In contrast to other existing models, our within-host SARS-CoV-2 model does not account for the time delay between a cell being infected and it becoming infectious, i.e., actively producing virus, termed eclipse phase (Figure 1A).Considering that the acute infection is on average 20 days [6], an eclipse period of about 6 hours [15] is short in comparison.Hence, we do not expect the omission of the eclipse phase to change the results qualitatively.Moreover, the immune response in the model is summarized by the dynamic variable B, which regulates the overall infectivity.There is a large variation on how modelling approaches incorporate the immune response and its effects on the viral dynamics of SARS-CoV-2 throughout literature [3,14,15,20].
Our model uses an approach, which properly captures more of the dynamic features of the immune response, such as waning immunity, while at the same time describing the short-term viral load dynamics as clinically observed (Figure 2A).Therefore, we believe our model assumptions to be justified.Furthermore, our aim was to develop a within-host SARS-CoV-2 model, which not only describes the acute infection phase, but also the post-acute infection dynamics and reinfection.We therefore assumed a simple optimization approach for describing the individual viral load dynamics and performed a qualitative comparison between our model fits and the original fits of Ke et al.
(Figure 2).Finally, it must be noted that the possibility of reinfection with the same virus variant is not accounted for in the model.While, in theory this is possible as a result of waning immunity, in reality the slow nature of waning combined with the rapid evolution of SARS-CoV-2 makes is unlikely for the variant causing the primary infection to be around by the time immunity is low enough to allow reinfection with it.
Although developed for the infection of cells in the upper respiratory tract the model can be easily adapted to account for the lower respiratory tract and extended to assess the effects of treatment options and drug resistance development.
Overall, the model presents an advance in the description of the SARS-CoV-2 infection dynamics within humans capturing the full dynamic infection process from initial infection to the clearance of the virus to reinfection and the corresponding short-term and persisting immune response across multiple infections.This modelling approach should be of interest for clinical use when quantitatively describing the within-host SARS-CoV-2 infection and developing treatment options seeking optimal treatment design.

Methods
Determining immune threshold B thres We defined immune threshold B thres as the minimal immune response required to counter a repeated spreading of the same virus variant within an individual.To determine B thres , we first derived the basic reproduction number R 0 , which signifies the number of secondary virus particles generated by the infection of a single virus when introduced into a fully susceptible cell population, hence, S = S 0 .If R 0 > 1, more virus particles are generated each generation leading to an increase of the viral population.If R 0 < 1, less virus particles are generated each generation letting the infection in the individual run out over time.To calculate R 0 , we first determined the Jacobian of the infected subsystem including compartments I and V : where T corresponds to the transmission and Σ to the transition matrix [5].The next-generation matrix NGM is then given by whose dominant eigenvalue is the reproduction number R, where R is dependent on B and hence time t (Figure 6): The reproduction number R can be split into , where (i) represents the number of viral particles produced during the life span of an infected cell and (ii) denotes the number of susceptible cells infected during the life span of a free viral particle.At the time of infection, the immune response is naive, such that the initial reproduction number R 0 without immune response (B = 0) is given by: The value B thres is then determined so that R = 1 for B = B thres .In the case p V ≤ d I , it is easy to show that R < 1 for any B ∈ [0, 1].This means that infection is not possible in this case, therefore we futher assume that p V > d I .Assuming the susceptible cells S reach their homeostatic level S 0 at the end of an infection, B thres in ODE system (1) is given by We assume B thres > 0, to ensure the possibility of infection within the individual.

Steady states
The steady states are defined as the states for which the right-hand sides of the ODE system (1) are all equal to zero.We directly see that the pre-infection steady state is given by (S * , I * , V * , B * ) = (S 0 , 0, 0, 0).
To identify the stability of the disease-free steady state, we first considered the Jacobian of the full ODE system described by (1), and determined the Jacobian at the disease-free steady state: We then identified the eigenvalues from with λ 1 = −d S < 0 and λ 2 = d B B thres > 0. From the remaining quadratic equation , where c * < 0, as < 1 according to the definition of B thres .Then the eigenvalues are , with a * = 1, and we get λ 3 < 0 and λ 4 > 0. Overall, this signifies that the disease-free steady state is unstable.Any initial amount of virus activates the immune response B, which then approaches B thres .Similarly, upon successful infection and recovery the post-infection steady state is given by (S * * , I * * , V * * , B * * ) = (S 0 , 0, 0, B thres ), where the immune response reaches memory B thres countering reinfection of the same virus variant.The Jacobian at (S * * , I * * , V * * , B * * ) is given by From this equation we get eigenvalues λ 1 = −d S < 0 and λ 2 = −d B B thres < 0. For eigenvalues λ 3,4 we need to solve the remaining quadratic equation When inserting B thres from (4), we see that c * * = 0 and hence, λ 3 = −(d V +β 0 (1−B thres )S 0 +d I ) < 0 and λ 4 = 0. Thus, the investigation of stability of the steady state X * * := (S 0 , 0, 0, B thres ) is not determined by its linearization and requires higher order analysis.According to the center manifold theorem (see e.g.Theorem 5.1 in [18]), there is an one-dimensional invariant manifold (a curve) C passing through X * * tangent to the eigenvector l (λ4=0) , corresponding to the eigenvalue λ 4 = 0.The stability of the whole system is determined by its stability on the curve C. On this curve, the system can be either be instable, semi-stable (on one part of the curve only) or stable.The eigenvector l (λ4=0) is given by The analytical description of the center manifold and the investigation of its stability is non-trivial.Hence, we rely on the numerical results, which clearly show local stability of system (1).In such a case, C is a so-called "slow manifold": the convergence to the steady state is not exponential but is at most so fast as 1/ √ t when time t increases to infinity.This is also consistent with the numerical results.There exists a third unique steady-state, however, not within the positive orthant, which is invariant.Since all eigenvalues are real and the trajectories of the system converge locally to the steady state X * * , the latter is globally stable in the positive orthant.Generally, this steady state can be approached asymptotically in one of two ways: either from above, when the short-term acute immune response is strong enough to drive B(t) > B thres in the initial phase of infection or from below.From the ODE system (1) we see that ∂B ∂t = 0 when (i) ( V , B) = (0, B thres ) For the latter to make sense biologically V ≥ 0, and hence it is required that B ≥ B thres .Whether the short-term acute immune response is strong enough to lead to B ≥ B thres depends on the viral load V and parameter p B .The definitions of B and B thres assure that the infection is fully and permanently cleared in the model, which is in line with current literature on SARS-CoV-2 for nonimmunocompromised individuals.Without immune threshold at B thres , the model would lead to a minimal but chronic infection with periodic spreading of the virus within the individual.

Data and pre-processing of single-individual SARS-CoV-2 infection dynamics
During fall 2020 and spring 2021, Ke et al. collected daily nasal samples for up to 14 days of all faculty, staff and students of the University of Illinois at Urbana-Champaign, who either (i) reported a positive quantitative reverse transcription polymerase chain reaction (RT-qPCR) result in the past 24 hours or (ii) were within five days of exposure to someone with a confirmed positive RT-qPCR result, while having tested negative for SARS-CoV-2 in the previous seven days [14].These criteria ensured a large-scale, high-frequency screening of early SARS-CoV-2 infection dynamics.In total, Ke et

Parameterization of the within-host SARS-CoV-2 model
To describe the individual SARS-CoV-2 infection dynamics, we parameterized the model as shown in Table 2. Initial values and the rate constants for five out of the eight parameters were taken from literature and assumed to be the same across individuals.The initial number of susceptible cells S 0 was derived by Ke et al. [14].Otherwise, we assumed the system to be initialized by a single successfully infected cell.Assuming that the system was in equilibrium during the initial disease-free state, we set the susceptible cell production rate constant p S to S 0 × d S , such that the in-and outflow of S is the same.The rate constant values for the infection rate β 0 and the death rate of an infected cell d I were taken from Ke et al. and are the mean rate constants of their estimated individuallevel rate constants [14].In the absence of an acute viral infection, the death rate constant of a susceptible cell d S is low, measured in the order of 0.02-0.03/day [28].However, this rate constant is estimated to increase during acute infection [25].For simplicity, we assumed a constant death rate constant between the lower and upper estimates.The value for the viral clearance rate constant d V was adopted from other studies [19,31].To account for the heterogeneity between individuals, constants of p V , p B , and d B were not within ±2 × standard deviations (stds) of their respective standardized estimated distributions.For these individuals, the CN nasal swab samples have either not been collected during early infection, such that infection dynamics before peak viral load are not well captured by the model, or where the CN values did not include sufficient information to deduce the dynamics of long-term immune response level formation.We ensured that each of the three remaining standardized distributions were standard normally distributed by performing a onesample Kolmogorov-Smirnov test employing the kstest function of MATLAB [12].The function kstest tests the null hypothesis that the parameter sample is drawn from a standard Gaussian distribution.As we required all three statistical tests to not reject the null hypothesis, we corrected for three-fold testing by applying the Bonferroni correction and adjusted the significance level of 0.05 to 0.05 3 = 0.0167 [1].The p-values for the standardized estimated distributions of p V , p B , and d B are 0.48, 0.99, and 0.19, respectively.By drawing from this standardized multivariate-Gaussian distribution, we maintained the correlation between the parameters as estimated previously and assured that their values were within a plausible range.In total, we drew 50 rate constant triplets out of the resulting standardized multivariate-Gaussian distribution and used these re-transformed rate constants to simulate 50 nasal viral infection dynamics for up to 20 days post infection (Figure 3B).

Comparison of key features
We compared five key features of the fitted and simulated infection dynamics, namely (i) the peak viral load, (ii) B thres , (iii) the number of days from infection to peak viral load, (iv) the number of days from peak viral load to undetectable viral load, and (v) the number of days from infection to B ≥ B thres (Figure 3B).For the fits and simulations which did not reach undetectable CN values or for which B ≱ B thres before 20 days post infection, we set these values to the maximum of 20 days.For each of the key features, we compared the resulting distributions from the fitted and simulated infection dynamics by performing a two-sample Kolmogorov-Smirnov test as provided by MATLAB function kstest2.The function kstest2 tests for the null hypothesis that both parameter samples come from the same continuous distribution.As we tested for a single universal null hypothesis, where all five alternative null-hypotheses were required to not be rejected, we corrected for fivefold testing by applying the Bonferroni correction and adjusting the significance level of 0.05 to 0.05 5 = 0.01 [1].

Reinfection
We differentiated between two different variants employing distinct modes of reinfection: (i) increased infectivity or (ii) immune escape.Increased infectivity was introduced into the model by an increased infection rate constant β0 > β 0 and B thres was modified to according to its definition (equation ( 4), Figure 5A).If β0 > β 0 , then Bthres > B thres .Similarly, immune escape was introduced into the model by an initially decreased level of the variant-specific immune response, B(0) < B( t), where B is the immune response of the new virus variant and t is the time of reinfection (Figure 5B).For this virus variant, the immune threshold remains the same, such that Bthres = B thres .For the simulations, we assumed individual #432192 to have been exposed to one of the variants, leading to a single successfully infected cell at 90 days post primary infection.
We then simulated the ODE system of (1) for the new virus variant for another 60 days.The initial values of the number of susceptible cells S and measured virus V , as well as the immune response B for the simulations of reinfection, were taken from the long-term fit of individual #432192 at 90 days since primary infection.We set the increased infection rate arbitrarily to β0 = 2β 0 and immune escape to 50%, such that B(0) = 0.5B( t).

Implementation and code availability
There is no original data underlying this work.Only previously published data was used for this study [14].The MATLAB code corresponding to this manuscript will be made available upon acceptance of the manuscript.The analysis was performed with MATLAB 2023a.

Figure 1 :
Figure 1: Schematic presentation of the within-host SARS-CoV-2 infection model and diagrammatic representation of the immune response during primary and re-infections.(A) Susceptible cells (S) are produced at rate p S and become infected at rate β 0 r(t), where β 0 is the infectivity rate, r(t) = 1 − B(t) is the immune response effect function and B(t) is the immune response.Infected cells (I) in turn produce free virus (V ) at rate p V .Susceptible and infected cells die at rates d S and d I , respectively, and virus is cleared at rate d V .Immune response B depends on the viral load V , is activated with rate p B , and decreases with rate d B .(B) Upon a primary SARS-CoV-2 infection, the acute short-term immune response detects, contains, and eliminates the virus.After viral clearance, immune protection against SARS-CoV-2 typically weakens (wanes) over time.In the case of reinfection with a different variant, the immune response is re-activated.For a different variant with immune escape (yellow), the effectiveness of the initial immune response against the different variant suffers a drop due to the variant's escape properties.For a different variant featuring higher infectivity (blue), the immune response rises to higher levels in order to successfully control the infection (blue).The dotted line represents the immune response without reinfection.

Figure 2 :
Figure 2: Model captures the individual-specific viral load dynamics of a primary SARS-CoV-2 infection.(A) Model fits of our model (orange line) and Ke model (gray line) to nasal CN values (black dots, measured by Ke et al.) of 12 randomly selected individuals.The dashed line represents the detection threshold of the RT-qPCR method used by Ke et al. to determine viral presence in the nasal swab sample and is set to CN = 42.Dots on the dashed line denote measurements below the detection threshold.(B) Comparison scatter plot of the mean squared errors (MSEs) of our model fits to fits from Ke et al. (MSE Ke ) for all 56 individuals.The dashed line indicates the diagonal representing equal MSEs between model fits.All points above the diagonal represent individuals for which the Ke model fits lead to smaller MSEs than our model (gray shaded area), while all points below the diagonal represent individuals, for which our model fits lead to smaller MSEs than the Ke model (orange shaded area).(C) Comparison of log 10 MSEs of model fits to all 56 individuals according to the Bland-Altman method.The mean of the log 10 MSE differences between our model fits and those from Ke et al. is -0.035 (black dashed line) and the mean of the log 10 MSE differences ±2 × standard deviations (stds) are -0.34 and 0.27, respectively (black dotted lines).All points above zero represent individuals for which the Ke model fits lead to smaller log 10 MSEs than our model (gray shaded area), while all points below zero represent individuals, for which our model fits lead to smaller log 10 MSEs than the Ke model (orange shaded area).

Figure 3 :
Figure 3: Model reproduces key features of SARS-CoV-2 infection dynamics.(A) Median, 10 th to 90 th percentiles of 42 fits of our model to nasal CN values (Ke et al., orange) and 50 simulated infection dynamics (purple) using our within-host model.The dashed line represents the detection threshold of the RT-qPCR method used by Ke et al. to determine viral presence in the nasal swab sample and is set to CN = 42.(B) Distributions of the five key features describing the underlying infection dynamics: (i) log 10 peak viral loads (VL), (ii) B thres -values and durations from (iii) infection to peak VL, (iv) peak VL to undetectable VL, defined as a CN value ¿ 42, and (v) infection to B > B thres across the fits (orange) and simulations (purple).The distributions of key features across fitted and simulated infection dynamics are statistically comparable with p-values 0.37, 0.51, 0.42, 0.52, and 0.45, respectively (two-sample Kolmogorov-Smirnov test correcting for five-fold testing according to the Bonferroni correction).The scatter is consistent between fits and simulations across sub-panels.The median values per feature are highlighted by the black bars.

Figure 4 :
Figure 4: SARS-CoV-2 infection dynamics following acute infection.SARS-CoV-2 infection dynamics for the populations of susceptible and infected cells S and I, respectively, free virus V , and the immune response B for all 56 fitted individuals following acute infection.The susceptible cells are normalized per individual by the initial number of susceptible cells S 0 , while the number of infected cells and free virus particles are normalized by their maximal values during infection.

Figure 5 :
Figure 5: Model fits and simulation of reinfection with virus variants with two different changes in properties: increased infectivity and immune escape.(A) Schematic presentation of the adapted model parametrization accounting for a different virus variant with increased infectivity rate constant β0 (here, β0 = 2β 0 ).Model fit of the viral load V and immune response B during and post-acute infection of a representative individual is shown in orange.The vertical dotted black line indicates the day of contact with the respective virus variant.The continued simulation demonstrates an increase in the viral load of the variant Ṽ , and hence, reinfection (blue line).The variant-specific immune response dynamics of B and corresponding variant-specific immune threshold Bthres are shown in the right panel (blue line and blue dashed line, respectively).The blue shaded area corresponds to the duration for which B < Bthres and indicates the period at which the new virus variant with increased infectivity can spread within the host.(B) Schematic presentation of the adapted model parameterization accounting for a different virus variant with immune escape reflected by a decreased initial variant-specific immune response B (here, B(0) = 0.5B(90)).Like in (A), the model fit of viral load V and immune response B during and post-acute infection of a representative individual is shown in orange.The vertical dotted black line indicates the day of contact with the respective virus variant.The continued simulation demonstrates reinfection (yellow line).The variant-specific immune response dynamics of B and corresponding variant-specific immune threshold Bthres are shown in the right panel (yellow line and dashed line, respectively).The yellow shaded area corresponds to the duration for which B < Bthres and indicates the period at which the new virus variant with immune escape can spread within the host.

Figure 6 :
Figure 6: Changes of the reproduction number R over the immune response B (A) and over the course of an infection (B).The parameter values are taken from Table 2, with p V = 200, p B = 10 −8 , and d B = 10 −2 .The horizontal dashed line denotes R = 1 and the vertical dotted line denotes the corresponding value of B, defined as the immune threshold B thres .
al. reported the cycle number (CN) values of nasal swab samples over time for 60 individuals.Due to very low or undetectable viral loads, we removed four out of the 60 individuals prior to the analysis, as also done by Ke et al.Time was reported relative to the day at which the maximal CN value was measured.The model is, however, initialized at infection.To compare measured and simulated CN values on the same time scale, we assumed viral load to peak six days post infection, as reported throughout literature for early SARS-CoV-2 variants [16, 26].Hence, we shifted the relative time scale of the CN values measured by Ke et al. by six days and removed all measurements prior to the assumed day of infection.Moreover, to directly compare CN values and simulated viral loads, we made use of the CN value-to-viral load calibration determined by Ke et al. and given by CN = − log 10 V − 11.35 0.25 .
[15]d S , β 0 , d I , p V , d V , p B , and d B , where all parameters are restricted to positive values for biological interpretability.V as the viral production rate times the sampled virus[15].To account for the heterogeneity between individuals, we estimated three of the eight model rate constants, namely the viral production rate constant p V , the activation rate constant of the immune response p B , and the waning rate constant d B , at an individual-specific level.All other rate constants were assumed to be fixed and the same across individuals.We observed a good fit of our model to the individual-specific CN values with similar dynamics as predicted by the model ofKe et al. (Figure

Table 1 :
Median values of key features of infection dynamics: output from our model (fitted and simulated) alongside corresponding values reported by earlier studies.