Development of a novel mathematical model that explains SARS-CoV-2 infection dynamics in Caco-2 cells

Experimental data

The main experimental data were obtained from the recent work of Shuai with co-authors (Shuai et al., 2022). Specifically, human colorectal adenocarcinoma cell lines Caco-2 (30,000 cells) were incubated 2 h with SARS-CoV-2 at multiplicity of infection (moi) of 0.1 ( $t=0$ stands for the start of incubation). After incubation cells were washed and media was replaced. Next, SARS-CoV-2 titers were assessed by TCID ${}_{50}$ assay at three points: $t=8,24,48$ h post infection (hpi). The number of subgenomic envelope viral RNA (sgE) relative to the human housekeeping gene GAPDH was measured by RT-qPCR at four time points ( $t=2,8,24,48$ hpi). Thus, RT-qPCR data included only intracellular viral RNA abundance and did not take into account the genomic RNA in virus particles. The data was obtained from the “Source Data Fig. 1” and “Source Data Extended Data Fig. 1” files deposited alongside with the manuscript (Shuai et al., 2022). The results of cell viability assays ( $t=2,24,72,120$ hpi) and the quantification of infected cells by N protein staining at $t=24$ hpi under the same experimental conditions were obtained from two other papers of the same group: “Fig. 3A” in Chu et al. (2020) and “Fig. 2E” in Shuai et al. (2020), respectively. The absence of cytopathic effect during SARS-CoV-2 infection in Caco-2 cells was also reported in multiple studies by other research groups (Wurtz et al., 2021; Zupin et al., 2022; Alexander et al., 2020; Bartolomeo et al., 2022; Saccon et al., 2021), including results for the Delta (Mautner et al., 2022) and the Omicron (Dighe et al., 2022) variants.

Standard ODE models

Usually, differential equation models of within-host viral dynamics of SARS-CoV-2 take into account the number of SARS-CoV-2 viral particles (V) and the numbers of uninfected (U) and infected (I) cells at the time moment $t$ after infection (Hernandez-Vargas & Velasco-Hernandez, 2020; Fadai et al., 2021; Du & Yuan, 2020):

(1) $$\{\begin{array}{l}{\displaystyle \frac{\text{d}U}{\text{d}t}=-\beta UV}\\ {\displaystyle \frac{\text{d}I}{\text{d}t}=\beta UV-\delta I}\\ {\displaystyle \frac{\text{d}V}{\text{d}t}=pI-cV}\end{array}$$

Here uninfected cells become infected proportional to the number of uninfected cells and virus particles at rate $\beta$, infected cells produce new viral particles at rate $p$, infected cells and viral particles decay at per capita rates $\delta$ and $c$, respectively (the latter implicitly includes immune system response). At the initial moment of time a healthy host ( $U(0)=U_0$, $I(0)=0$) becomes infected with $V(0)=V_0$ viral particles. A similar model was used by Bernhauerová et al. (2021) to explain in vitro SARS-CoV-2 infection data. However, we showed that the solutions of these differential equations contradict to the cell viability data (see Results).

The models proposed

We used two models of the SARS-CoV-2 infection process: the “discrete” model and the “continuous” model. The discrete model in essence implements the approach proposed by Doob et al. and Gillespie (see e.g. Gillespie, 1977), but is additionally endowed with a number of continuous parameters, i.e. the ones that take real values, so possibly a more correct term is “semi-discrete”. The continuous model can be viewed as the limit of the discrete model; it is the system of integro-differential equations in the variables $U_0$ (“regular” healthy cells), $U_1$ (resistant cells), $Cyt$ (cytokines), I (infected cells), V (free virions) and R (intracellular viral RNA).

Continuous model

The continuous model can be viewed as the limit of the discrete model described above. The model is the system of integro-differential equations in the variables $U_0$ (“regular” healthy cells), $U_1$ (resistant cells), $Cyt$ (cytokines), I (infected cells), V (free virions) and R (intracellular viral RNA):

$$\begin{array}{rl}& \{\begin{array}{lll}{\displaystyle U_0(t{)}^{\mathrm{\prime }}}& =& {\displaystyle -\beta U_0(t)V(t)-{\beta }_{cyt}Cyt(t)U_0(t)}\\ {\displaystyle U_1(t{)}^{\mathrm{\prime }}}& =& {\displaystyle -\frac{\beta }{k}{U}_1(t)V(t)+{\beta }_{cyt}Cyt(t)U_0(t)}\\ {\displaystyle Cyt(t{)}^{\mathrm{\prime }}}& =& {\displaystyle p_{c}I(t-{\mathrm{\Delta }}_{tcyt})-{\beta }_{cyt}Cyt(t)U_0(t)}\\ {\displaystyle I(t{)}^{\mathrm{\prime }}}& =& {\displaystyle \beta (U_0(t)+U_1(t)/k)V(t)}\\ {\displaystyle V(t{)}^{\mathrm{\prime }}}& =& {\displaystyle {\int }_0^t(\beta (U_0(t-x)+(U_1(t-x))/k)V(t-x))p(x)dx-}\\ & & {\displaystyle -\frac{\beta }{k}(I(t)+U_1(t))V(t)-\beta U_0(t)V(t)}\\ {\displaystyle R(t)}& =& {\displaystyle {\int }_0^t(\beta (U_0(t-x)+U_1(t-x)/k)V(t-x))\sigma (x)dx}\end{array}\end{array}$$

Here $p(t)$ is the first derivative of the number of virions produced by an infected cell (evaluation of this function is described in Appendix S2), $\sigma (t)$ is the total concentration of viral RNA produced by a cell up to the time $t$ after infection ( $\sigma (t)=2^{t/{\mathrm{\Delta }}_{tRNAdouble}}$). The first equation means that “regular” healthy cells become infected with intensity proportional to the product of the number of healthy cells and the number of free virions or become resistant with intensity proportional to the product of the number of healthy cells and the number of cytokines. The second equation means that resistant cells are produced from “regular” healthy cells and can be infected with intensity proportional to the product of the number of resistant cells and the number of free virions (note that the coefficient in the second equation differs from the coefficient in the first equation; in terms of Fig. 2 this coefficient is ${\beta }^{\mathrm{\prime }}=\beta /k$). The third equation is devoted to cytokines: new cytokines are produced by infected cells with a certain delay, and existing cytokines are used to make healthy cells resistant. The fourth equation describes infected cells which can be produced either from “regular” healthy cells or from resistant cells. The fifth equation is devoted to release of free virions. The first summand in essence means that infected cells produce free virions with time-dependent intensity $p(\tau )$ (note that the first factor under the integral sign is the right-hand side of the fourth equation). Speaking in more detail, this summand can be rewritten in the form

$$\begin{array}{rl}& {\displaystyle {\int }_0^t(\beta (U_0(t-x)+(U_1(t-x))/k)V(t-x))p(x)dx={\int }_0^t{I}^{\mathrm{\prime }}(t-x)p(x)dx;}\end{array}$$free virions are produced by cells infected at different time points. Every cell infected in the interval $[t-x,t-x+dt]$ (the number of such cells is $I^{\mathrm{\prime }}(t-x)dt$) at the interval $[t,t+dt]$ emits $p(x)dt$ virus particles. Thus the total number of virions produced by cells infected in the interval $[t-x,t-x+dt]$ is $I^{\mathrm{\prime }}(t-x)dtp(x)dt$. Hence, it only remains to integrate these values. The second summand shows that free virions are used to infect healthy cells but also may enter cells that are already infected. Finally the sixth equation computes concentration of viral RNA inside infected cells similarly to production of free virions in the fifth equation. The sixth equation was added to the system solely for the purpose of fitting to the existing data, as the variable $R(t)$ is not used in other equations.

Note that the sum of the first, the second and the fourth equation has zero in the right-hand side, thus the total number of cells is constant, i.e. cells do not die.

Computational approach for solving the proposed system of integro-differential equations is presented in Appendix S2.

Selection of parameters of the proposed models

Our models operate with four types of parameters. The values of parameters of the type “Fixed” are extracted from literature or are directly inferred from data. The values of parameters of the type “Global” are assumed to be equal between the models corresponding to different virus types and are tuned using the whole dataset. One part of the “Global” parameters ( ${\beta }_{cyt}$, $k$) is exclusively related to the properties of the target cells. The other “Global” parameters ( $n_{l2V}$, $p_{Vir}$, ${\mathrm{\Delta }}_{tRNAdouble}$, ${\mathrm{\Delta }}_{tlatent}$) are shared by the WT and the Delta viruses because we assume that the corresponding biological values could not be significantly altered by such a low number of mutations in viral genomic RNA. The values of parameters of the type “Local” are tuned using data corresponding to a single virus type. Finally, parameters corresponding to the implementation of the models are of the type “Internal”.

Data description and motivation

We started with the qualitative analysis of previously published dataset of Caco-2 cells infected with the Wuhan (WT) and the Delta variants of SARS-CoV-2. Specifically, the measurements of infectious viral particles in medium, viral subgenomic RNA in infected cells and cell viability data were available for several time points after infection (Fig. 1). We also used the data on the percentage of infected cells 24 h post infection (see Materials and Methods for the details). Based on these data we made three key observations.

1. Caco-2 cells (including infected ones) stay viable throughout the whole experiment.

2. Both virus titers and intracellular viral RNA growth significantly slow down over time: 5 folds difference in the WT virus titers growth rates between 8–24 and 24–48 hpi intervals, 4.7 folds in the case of the WT viral RNA growth, 2.9 folds for the Delta virus titers, and 5.7 folds for the Delta viral RNA.

3. The fraction of infected cells after 24 h post infection is relatively small (approximately 6%).

Since the number of uninfected cells is sufficiently large throughout the whole experiment (point 3), exhaustion of healthy susceptible cells cannot explain the decrease in infection rate (point 2). This observation along with the absence of cytopathicity (point 1) led us to the hypothesis on the exhaustion of some key cellular resources needed for virus particle assembly and release (e.g., lipids). We also hypothesized that the primary reason for the lack of infection in the majority of cells is connected to cell immunization (innate immune response).

The existing standard mathematical models (see Materials and Methods) were not applicable for explaining such data for several reasons. First, these models operate with the numbers of cells (uninfected, infected) and viral particles. Thus, measurements of intracellular viral RNA could be used for parameter inference only if we assume proportionality of viral RNA to the number of infected cells (which, as we believe, is a too strong assumption). Alternatively, one can infer parameters using only infectious virus titers data which ultimately results in unidentifiability of some parameters (Bernhauerová et al., 2021). Second, the absence of a cytopathic effect results in solution properties which contradict experimental data. Namely, in the case of the model with three variables (U, I, V) all cells become infected, which implies an unbounded increase in the number of viral particles, since they are produced by infected cells at a constant rate (see S4 Appendix for an accurate proof).

A new mathematical model of SARS-CoV-2 infection in vitro

Given the aforementioned observations, we developed a mathematical model which included a possibility of cell exhaustion, direct modeling of intracellular viral RNA and the innate immune response. The model is schematically illustrated in Fig. 2. A healthy cell gets infected with rate constant $\beta$ upon contact with a productive virion. An infected cell produces virions after latent period ${\mathrm{\Delta }}_{tlatent}$ proportional to the number of viral genomes and cellular resources. Once a cell becomes exhausted (i.e., resources for virion assembly/release run out), it can no longer produce new viral particles. The described scheme was initially explicitly implemented in the discrete model, while for the continuous differential equation model the virus production rate was modeled by the associated function $p(t)$ (Appendix, Fig. S1; see Materials and Methods for details). An infected cell also produces cytokines with the rate constant $p_c$ with a delay ${\mathrm{\Delta }}_{tcyt}=\tau$. Cytokines then enter healthy cells with rate constant ${\beta }_c$ and make them more resistant to the virus; infection rate constant for resistant cells is ${\beta }^{\mathrm{\prime }}<\beta$.

Parameter fitting and comparative analysis

For both WT and Delta viruses we optimized parameters of the mathematical model so that solutions fit the experimental data (virus titers, subgenomic RNA, percentage of the infected cells at $t=24$ hpi). To simplify the computational procedure and improve interpretability, we assumed that most of the model parameters were the same for two SARS-CoV-2 variants. The exceptions were virus cell entry rate $\beta$ and cytokine production intensity $p_c$. For both variants the solutions explained the data well (Fig. 3; Appendix, Fig. S2). In both cases only a fraction of infected cells were active (i.e., producing new viral particles) at $t=48$ hpi. Notably, the discrete model allowed us to explicitly quantify the number of active infected cells (Fig. 3), which was not a case in differential equation model. Solutions for the remaining variables were near-identical between both models (Appendix, Fig. S2).

Note that graphs for the number of virions contain flat sections. The reason for the first flat section is the latent period (approximately $7$ h, see the parameter ${\mathrm{\Delta }}_{tlatent}$ in Table 1): all virus particles are removed by washing, and new particles have not been emitted yet. The following flat segments are the reflections of the first one, since virus production period of an infected cell is short due to resource exhaustion (in terms of the continuous model the function $p(x)$ in the fifth equation of the system becomes equal to $0$; see Appendix, Fig. S1), and new infected cells are in the latent period.

The key difference between replication profiles of two viruses was a significant reduction in the number of intracellular viral subgenomic RNA in the Delta-infected cells compared to the WT: 2.41 folds at 24 hpi, 4.2 folds at 48 hpi (Mann–Whitney’s U-test $p=2.16\times {10}^{-3}$ for both time points). A similar trend was also observed for infectious virus titers data, though the difference at 24 hpi (7.8 folds, $p=0.029$) was higher than at 48 hpi (3.1 folds, p = 0.057). Our model captured these data well, extrapolating the same trends to the absolute (non-observable) numbers of viral particles and infected cells (Fig. 3). Both cell entry and cytokine production rates had different values as a result of parameter fitting to the WT and Delta data. Specifically, the rate of virus entry into Caco-2 cells was 1.7-fold higher for the WT compared to the Delta variant. Conversely, cytokine production rate by an infected cell was higher for the Delta variant (1.4-fold). Confidence intervals for both parameters were not overlapping.

Full list of parameter values is presented in Table 1. We also verified that if all parameters except one are equal for WT and Delta, then approximation quality essentially drops (the least possible value of $Err$ in this case is 1.8 times greater than for the case considered above).