Intracellular Life Cycle Kinetics of SARS-CoV-2 Predicted Using Mathematical Modelling

Abstract

SARS-CoV-2 infection represents a global threat to human health. Various approaches were employed to reveal the pathogenetic mechanisms of COVID-19. Mathematical and computational modelling is a powerful tool to describe and analyze the infection dynamics in relation to a plethora of processes contributing to the observed disease phenotypes. In our study here, we formulate and calibrate a deterministic model of the SARS-CoV-2 life cycle. It provides a kinetic description of the major replication stages of SARS-CoV-2. Sensitivity analysis of the net viral progeny with respect to model parameters enables the identification of the life cycle stages that have the strongest impact on viral replication. These three most influential parameters are (i) degradation rate of positive sense vRNAs in cytoplasm (negative effect), (ii) threshold number of non-structural proteins enhancing vRNA transcription (negative effect), and (iii) translation rate of non-structural proteins (positive effect). The results of our analysis could be used for guiding the search for antiviral drug targets to combat SARS-CoV-2 infection.

1. Introduction

Human infection with SARS-CoV-2 presents a tremendous health problem. The within-host infection characteristics are characterized by extreme variability of the disease course ranging from asymptomatic infections to severe forms of COVID-19 with lethal outcomes [1]. This spectrum of pathogenicity is a result of the interaction of numerous processes and factors on multiple levels of realization: the single cell, tissues and organs, and the organism’s physiology [2]. The SARS-CoV-2 infection dynamics in a human organism is determined by the kinetics of infection of cells which express viral receptors (e.g., ACE2), by the activation of the intracellular defense, and by systemic immunophysiological reactions. The corresponding “virus–human organism” is a multiscale, multicomponent dynamical system. A comprehensive study of such a system requires the development of mathematical models that integrate the underlying biophysical, biochemical, and physiological processes [3]. Although a few mathematical models have been recently proposed to describe the within-host kinetics of SARS-CoV-2 infection [4,5,6,7,8,9], the degree of a mechanistic process resolution remains to be greatly enhanced.

From a scientific point of view, the intracellular replication of the viruses is realized as an intertwined set of biochemical reactions and biophysical transport processes. Being summarized in a schematic form, the theoretical abstraction of the system provides a conceptual platform for deriving a mathematical model. The mathematical descriptions can further be specified on the basis of the law of mass action, enzyme kinetics, diffusion laws, etc, following deterministic, stochastic, or hybrid approaches. The derived mathematical models which bear a direct resemblance to and coordination with the underling chemical and physical processes are considered as mechanistic models versus empirical (e.g., statistical) or black-box (e.g., neural-networks) type mathematical models. The mechanistic models provide a descriptive and analytical tool for studying the virus life cycle regulation and predicting its response to structural or parametric perturbations. In virology, the mechanistic models have been successfully developed and applied to study, in a quantitative way, the intracellular replication of HIV-1 [10], hepatitis B virus [11], influenza A virus [12], hepatitis C virus [13], and poliovirus [14].

The dynamic view on the balance between virus expansion and development of antiviral immune responses suggests that the kinetics of virus growth is one of the major determinants of the infection trajectory [15]. Existing mathematical models of SARS-CoV-2 infection consider the virus growth at the cell population level, ignoring details of the virus life cycle. The infection of permissive cells by the virus results in a number of dramatic changes of the infected cell physiology, such as modification of host protein synthesis [16], inhibition of the innate immune responses [17], and induction of programmed cell death [18]. In this study, we formulate a mathematical model of SARS-CoV-2 replication in infected cells. The development of the model follows a deterministic approach, similar to [10], in order to focus on the calibration of model parameters to reproduce some reference kinetics of the virus life cycle steps. By applying a sensitivity analysis of the model, we rank the model parameters according to their impact on the net secretion of virions by the infected cell. The most influential parameters can be considered as prospective targets for available antiviral drugs or novel treatments.

2. Methods

2.1. Model Development and Calibration

The kinetics of the corresponding biochemical reactions is described in the deterministic mathematical model introduced in Section 3. The system of ordinary differential equations (ODEs) is formulated using the law of mass action and Michaelis–Menten parameterization. Note that the deterministic version provides an initial step toward a stochastic description of virus replication.

We use a Michaelis–Menten-type description for the nucleocapsid formation and viral assembly processes to represent the following aspects of virus replication: (i) purified coronavirus virions containing mainly full-genome length viral RNAs, i.e., there is no saturation in the kinetics with respect to viral RNA [18], (ii) structural proteins (N and M) required for virion morphogenesis and recruitment of the virion structural components to the assembly site [1] scale this rate so that enough proteins must be present above a certain level to reach half the maximal rate of the ribonucleoprotein formation and virion assembly, and (iii) the least abundant structural proteins represent the kinetic bottleneck in the overall progeny production.

A major step in the mathematical model formulation is the estimation of parameters appearing on the right-hand sides of the model equations. The maximum likelihood approach represents a general framework to address the parameter estimation problem. However, the sparsity of the available kinetic data on the intracellular life cycle of SARS-CoV-2 motivates the implementation of the parameter estimation procedure, which we refer to as “model calibration”. It stands for an iterative process of parameter guessing and constraining the choice progressively by requiring consistency with all available data. In this study, these data mainly consist of previously estimated numerical values and/or numerical ranges of the molecular species variables at key time points, such as the beginning and end of the infection cycle, or of intermediate steps, as summarized in Table 1.

viruses-13-01735-t001_Table 1 Table 1

Time-dependent variables of the mathematical model characterizing the SARS-CoV-2 life cycle.

Variable	Meaning	Quantitative Characteristics
$\left[V_{free}\right]$	number of free virions outside the cell membrane	10
$\left[V_{bound}\right]$	number of virions bound to ACE2 and activated by TMPRSS2	1–10
$\left[V_{endosome}\right]$	number of virions in endosomes	1–10
$\left[gRN{A}_{(+)}\right]$	single strand positive sense genomic RNA	1–5
$\left[NSP\right]$	population of non-structural proteins	−
$\left[gRN{A}_{(-)}\right]$	negative sense genomic and subgenomic RNAs	10
$\left[gRNA\right]$	positive sense genomic and subgenomic RNAs	10,000
$\left[SP\right]$	total number of structural proteins $S+M+E$ per virion	$2000\in (1125,2230)$ [19,20,21]
$\left[N\right]$	N proteins per virion $\left[N\right]$	456 [21]; $1465\in (730,2200)$ [19]
$\left[\mathit{N-gRNA}\right]$	ribonucleocapsid molecules	−
$\left[V_{assembled}\right]$	assembled virions in endosomes	−
$\left[V_{released}\right]$	virus burst size	10–10,000 virions in 7 to 24 h [2,22,23]

2.2. Model Validation

To validate the calibrated model, we compared some predicted quantities with available experimental data, namely, (i) the kinetics of positive- to negative-sense vRNA ratio in MHV-infected cells [24], and (ii) the kinetics of SARS-CoV-2 replication in the cell cultures of Vero E6 cells [23,25]. As we cannot directly infer the absolute values of the number of released virions by an infected cell just because there are no such single-cell experiments, we used the relative measure of fold changes of titer levels in cell culture from the start of their exponential growth. This data was compared with the fold changes of $\left[V_{released}\right]\left(t\right)$ starting from the moment $t_s$ at which $\left[V_{released}\right]\left(t_s\right)=1$.

The experimental data presented in [24] characterize the kinetics of virus transcription and translation of cells infected with Murine coronavirus (MHV-A59). In particular, the ribosome profiling technique was used to generate high precision data on the production of the positive and negative-sense genomic and subgenomic viral RNAs. The overall kinetics of SARS-CoV-2 in Vero E6 cells is studied in [23]. To this end, one-step growth (MOI = 5) of three recombinant (i.e., icSARS-CoV-Urbani, icSARS-CoV-GFP, icSARS-CoV-nLuc) and one clinical strain WA1 was followed. The data on the kinetics of virus titers in the supernatants of cell cultures (PFU/mL) show that the growth phase until reaching the plateau takes about 12 to 19 h and is similar to all the above CoV variants. A similar study of SARS-CoV-2 replication features such as growth kinetics, virus titers, analysis of transcription, and translation in Vero E6 cells is presented in [25]. The virus isolate SARS-CoV-2 Australia/VIC01/2020 was used for infection of cells at MOI = 3. Intracellular viral RNA, protein synthesis and release of infectious viral progeny (by plaque assay) were quantified. The initial guess for the mathematical model parameters was specified using data from a broad spectrum of publications covering the structural and genetic properties of the SARS-CoV-2, the coronavirus transcription and translation, and the turnover of proteins and RNA in cells with the relevant references presented in Section 3.

2.3. Parameter Uncertainty Analysis

Through the process of model calibration, the point estimates of parameter values were determined, which were used in main simulations and are discussed in Section 3. For most of the model parameters, we were able to estimate the ranges of biologically plausible values based on literature data. For other parameters, we derived their ranges applying the parameter uncertainty analysis. To this end, we iteratively adjusted their parameter ranges so that they would include their point estimates while the uncertainty of the model output (progeny release kinetics) would be confined within the range from 10 to 10,000 virions, as specified in Table 1. For some parameters with wide literature-based ranges, we further narrowed them down to restrict the output uncertainty. These three categories of parameter ranges (based on literature analysis, based on uncertainty analysis, based on both types of analysis), as well as parameters which are fixed ad hoc, are reviewed in Section 3 accordingly. To quantify the output uncertainty, we employed the Latin hypercube sampling method to randomly sample $n=10$,000 combinations of parameters from their respective ranges. Then, the median and 5–95% confidence regions of the obtained ensemble of $\left[V_{released}\right]\left(t\right)$ trajectories were computed.

2.4. Sensitivity Analysis

Let $\mathbf{u}(t,\mathbf{p})$ be the model solution of the ODE initial value problem $\frac{d\mathbf{u}}{dt}=\mathbf{f}(t,\mathbf{u},\mathbf{p})$, $\mathbf{u}\left(0\right)={\mathbf{u}}_0$. In this study, we apply the methods of local sensitivity analysis to identify the parameters (and respective biochemical processes) that have the most strong impact on the characteristics of interest $\mathrm{\Phi }\left(\mathbf{u}\right(\mathbf{p}\left)\right)$, e.g., the cumulative number of released virions

${\mathrm{\Phi }}_{progeny}={\int }_0^T{k}_{release}\left[V_{assembled}\right]dt,$

the area under the curve for $\left[V_{released}\right]\left(t\right)$

${\mathrm{\Phi }}_{progeny}^{AUC}\left(\mathbf{p}\right)={\int }_0^T\left[V_{released}\right]dt,$

or the distance between the model output $g\left(\mathbf{u}\right(\mathbf{p}\left)\right)$ and the experimental data $g_{exp}\left(t_i\right)$

$\mathrm{\Phi }=\sum _i \left|\right|g\left(\mathbf{u}\left(t_i\right)\right)-g_{exp}\left(t_i\right){\left|\right|}^2.$

The sensitivity index of model parameter p is defined as $s_p=\frac{d\mathrm{\Phi }\left(p\right)}{dp}$. It indicates the influence of parameter variation on the output value $\mathrm{\Phi }$. To compare sensitivities and rank parameters by their impact, the sensitivity indices are usually normalized by parameter values: ${\widehat{s}}_p=p\frac{d\mathrm{\Phi }\left(p\right)}{dp}$. The sensitivity indices can be computed by solving the forward sensitivity ODE system or by solving the adjoint problem [26]. We follow the adjoint-based approach as previously described [10,26].

2.5. Software

The following packages in Julia language (https://julialang.org, accessed on 12 May 2021) were used to simulate and analyze the model: DifferentialEquations.jl (accessed on 12 May 2021) for numerical solution of the model, QuasiMonteCarlo.jl (accessed on 12 May 2021) for parameter uncertainty analysis, DiffEqSensitivity.jl (accessed on 12 May 2021) for local sensitivity analysis, PyPlot.jl (accessed on 12 May 2021) for visualizations. The scripts used to simulate and analyze the model are provided in the Supplementary Materials.

3. Results

We follow the latest view of the SARS-CoV-2 life cycle summarized in [1,17,18,27]. The key steps include: (i) cell entry, (ii) genome transcription and replication, (iii) translation of structural and accessory proteins, and (iv) assembly and release of virions (Figure 1). The set of time-dependent molecular species described in the model is listed in Table 1 with some quantitative characteristics of their abundance.

3.1. Mathematical Model of Intracellular SARS-CoV-2 Replication

3.1.1. Cell Entry

The entry stage of SARS-CoV-2 is split into three steps as represented in Figure 1:

• binding of the receptor-binding domain (RBD) of the viral S protein to the ACE2 receptor,

• priming by host cell surface protease TMPRSS2,

• fusion at the cellular or endosomal membrane followed by release and uncoating of the viral genomic RNA.

Binding of the virion to the cellular transmembrane protein ACE2, and entry and release of the viral RNA into the host cell are described by equations specifying the rates of changes of free-, receptor-bound, and fused virions, as well as the viral RNA genome in the cytoplasm:

(1)$\begin{array}{cc}\hfill \frac{d\left[V_{free}\right]}{dt}& =-k_{bind}\left[V_{free}\right]-d_V\left[V_{free}\right]+k_{diss}\left[V_{bound}\right]\hfill \end{array}$

(2)$\begin{array}{cc}\hfill \frac{d\left[V_{bound}\right]}{dt}& =k_{bind}\left[V_{free}\right]-(k_{fuse}+k_{diss}+d_V)\left[V_{bound}\right]\hfill \end{array}$

(3)$\begin{array}{cc}\hfill \frac{d\left[V_{endosome}\right]}{dt}& =k_{fuse}\left[V_{bound}\right]-(k_{uncoat}+d_{endosome})\left[V_{endosome}\right]\hfill \end{array}$

(4)$\begin{array}{cc}\hfill \frac{d\left[gRN{A}_{(+)}\right]}{dt}& =k_{uncoat}\left[V_{endosome}\right]-d_{gRNA}\left[gRN{A}_{(+)}\right].\hfill \end{array}$

Here, $\left[V_{free}\right]$ is the number of free virions outside the cell membrane, $\left[V_{bound}\right]$ is the number of virions bound to ACE2 and activated by TMPRSS2, $\left[V_{endosome}\right]$ is the number of virions in endosomes, and $\left[gRN{A}_{(+)}\right]$ is the number of ss-positive sense genomic RNA. The respective parameters of the above equations are described in Table 2. In estimating the degradation rate of virions in endosomes, we followed the assumptions presented for modelling IAV infection [12], i.e., about $50\%$ fail to release the viral genome. This gives $d_{endosome}=0.06$ h${}^{-1}$.

3.1.2. Genome Transcription and Replication

The SARS-CoV-2 virion consists of about a 30 kb strand of positive sense RNA coated with N protein and covered by a lipid bilayer containing spike S, membrane M, and envelope E proteins [18]. The model was calibrated to reproduce (i) the scale of viral proteins production corresponding to about 10 to 10,000 infectious virions per cell, (ii) the observed ratio of positive and negative sense viral RNA (genomic and subgenomic) [24], and (iii) the known ranges of the parameters of mRNA and protein turnover [28].

The released genomic RNA undergoes translation into viral polyproteins (pp1a, pp1ab) which generate via proteolysis 16 non-structural proteins (nsp1–16). They are operating to form the viral replication and transcription complex. In particular, a key step is the formation of nsp12, which encodes the RNA-dependent RNA polymerase (RdRp). The primary function of the RdRp replication complex is to generate a negative sense full-length genome and subgenomic RNAs. It has been established for MHV virus that synthesis of negative-sense RNA starts about 60 to 90 min post-infection and reaches a maximum at about 5 to 6 h [24,29]. The resulting set of negative-sense RNAs and the full-length antisense genome are working as templates for the synthesis of positive-sense genomic and subgenomic RNAs as shown in Figure 1. The total number of produced positive-sense viral genomes and subgenomic RNAs exceeds the number of negative-sense RNAs by 100 to 1000 fold [24,29].

We describe the abundance of the populations of non-structural proteins $\left[NSP\right]$, the set of negative sense genomic and subgenomic $\left[gRN{A}_{(-)}\right]$, and the set of positive sense genomic and subgenomic $\left[gRNA\right]$ with the following differential equations:

(5)$\begin{array}{cc}\hfill \frac{d\left[NSP\right]}{dt}& =k_{transl}{f}_{ORF1}\left[gRN{A}_{(+)}\right]-d_{NSP}\left[NSP\right]\hfill \end{array}$

(6)$\begin{array}{cc}\hfill \frac{d\left[gRN{A}_{(-)}\right]}{dt}& =k_{t{r}_{(-)}}\left[gRN{A}_{(+)}\right]{\theta }_{RdRp}-d_{gRN{A}_{(-)}}\left[gRN{A}_{(-)}\right]\hfill \end{array}$

(7)$\begin{array}{cc}\hfill \frac{d\left[gRNA\right]}{dt}& =k_{t{r}_{(+)}}\left[gRN{A}_{(-)}\right]{\theta }_{RdRp}-(k_{complex}{\theta }_{complex}+d_{gRNA})\left[gRNA\right]\hfill \end{array}$

where

(8)${\theta }_{RdRp}=\frac{\left[NSP\right]}{\left[NSP\right]+K_{NSP}},{\theta }_{complex}=\frac{\left[N\right]}{\left[N\right]+K_N}$

It is taken into account in Equation (5) that the non-structural proteins are translated only from the released genomic RNA. Similarly, the transcription of the negative-sense genomic and subgenomic RNA described by Equation (6) is determined by the original positive-sense genomic RNA.

The translation rate $k_{transl}=45$,360 nt/mRNA h${}^{-1}$ has been estimated for coronaviruses in [24]. The length of the RNA genome coding for $\left[NSP\right]$ proteins is about 21,000 nucleotides [30], hence $f_{ORF1}=1/21$,000. These estimates are consistent with the general range of protein synthesis rates: $(1,{10}^4)$ molecules/mRNA h${}^{-1}$ [28].

The degradation rate of NSPs $d_{NSP}=0.069$ h${}^{-1}$ is estimated as the geometric mean of the half-lives of proteins in cells. This result of our calibration procedure is consistent with the finding that the lifetime of most proteins is just a few hours [28].

The transcription rates of negative sense genomic and subgenomic RNAs is estimated to be $k_{t{r}_{(-)}}=3$ copies/mRNA h${}^{-1}$. This is consistent with the transcription rate of mRNAs which ranges from 1 to 100 copies per hour according to [28].

We assume the threshold for half-maximal rate of RdRp activity to be $K_{NSP}=100$ copies, taking into account that a small number of non-structural proteins is sufficient for enhancing the transcription of vRNAs.

It is known that an average half-life of mRNAs in cells of vertebrates is about 3 h [31] and ranges from 1 to 10 h [28]. Hence, the decay rate of mRNAs can be in between $[0.069,0.69]$ h${}^{-1}$. We use here the following value $d_{gRNA}=0.2$ h${}^{-1}$. In addition, we assume a smaller value for the degradation rate of negative sense vRNAs in double-membrane vesicles, $d_{gRN{A}_{(-)}}=d_{gRNA}/2=0.1$ h${}^{-1}$.

Quantitative analysis of RNA polymerase elongation provides the estimate of the transcription rate to be $\mathrm{46,080}\pm \mathrm{17,640}$ nt/h [32]. Thus, the basal rate of transcription is around $\mathrm{46,080}/\mathrm{30,000}$ nt/RNA h${}^{-1}$. In infected cells, however, the overall rate of viral RNA transcription is amplified (around a 1000-fold increase) while transcription of host RNAs is largely silenced [30]. Therefore, we estimate $k_{t{r}_{(+)}}=1000$ copies/mRNA/h, which matches the observed ratio of positive- to negative-sense viral RNAs [24,29] (see the model validation results at the end of this Section).

The kinetics of the nucleocapsid formation (i.e., viral RNA genome coated with N protein) resulting from the binding of N proteins and gRNA is characterized by the rate constant $k_{complex}$. This can be estimated from the binding data presented in [33]. The data indicate a fast kinetics with a characteristic time of ≈20 s. Hence $k_{complex}\approx 0.4$ h${}^{-1}$, taking into account that 1 virion consists of 38 ribonucleoprotein complexes each having about 12 N proteins [21], $n_N=38\times 12=456$.

The kinetics of formation of the nucleocapsid condensates from N proteins and gRNA has been studied in [34,35]. The minimum concentration of N proteins necessary to form condensates has been estimated to be about $3.3\mathsf{\mu }$M, and the concentration at which the formation slows down to be about $11\mathsf{\mu }$M. These concentrations correspond to approximately $1.5$ and 5 million molecules, if we take 768 fL as the mean cell volume of type II pneumocyte [2]. The estimate $K_N=5\times {10}^6$ molecules implies that the saturation takes place at the number of N proteins necessary to assemble around 11,000 virions ($K_N\approx n_N·\mathrm{11,000}$).

3.1.3. Translation of Structural and Accessory Proteins

The structural proteins S, envelope E, and membrane M are translated from the positive sense subgenomic RNAs at the endoplasmic reticulum (ER). They are described in the model by their total abundance $\left[SP\right]$. They are considered to interact together and assist in forming virus-like particles and budding of new virions from the ER and Golgi compartments (ERIGC) [17]. The structural nucleocapsid protein N is translated from subgenomic RNAs by cytosolic ribosomes and can enhance the formation of virus-like particles [27]. Its key function is to create nucleocapsid $\left[\mathit{N-g}RNA\right]$ by coating genomic RNAs. The number of N proteins per virion $\left[N\right]$ in coronaviruses is estimated to range from 730 to 2200 [19]. For SARS-CoV-2, however, the estimated number of N proteins per virion is $n_N=38\times 12=456$ [21].

The translation rates of $\left[N\right]$ and $\left[SP\right]$ proteins are described by the following two equations:

(9)$\begin{array}{cc}\hfill \frac{d\left[N\right]}{dt}& =k_{transl}{f}_N\left[gRNA\right]-k_{complex}{n}_N{\theta }_{complex}\left[gRNA\right]-d_N\left[N\right]\hfill \end{array}$

(10)$\begin{array}{cc}\hfill \frac{d\left[SP\right]}{dt}& =k_{transl}{f}_{SP}\left[gRNA\right]-k_{assemb}{n}_{SP}{\theta }_{assemb}\left[\mathit{N-gRNA}\right]-d_{SP}\left[SP\right]\hfill \end{array}$

where

(11)${\theta }_{assemb}=\frac{\left[SP\right]}{\left[SP\right]+K_{{\mathrm{V}}_{\mathrm{rel}}}{n}_{SP}}$

The scaling parameter $f_N=1/1200$ accounts for the length of the N-coding RNA which is about 1200 nucleotides (400 aa) [36]. Likewise, $f_{SP}=1/10$,000, as the estimated RNA length for the structural proteins S, E, and M is about 10,000 nucleotides.

The degradation rate of N protein $d_N=0.023$ h${}^{-1}$ is estimated using the database [37]. The degradation rate of the mixture of other structural proteins is evaluated to be around $d_{SP}=0.044$ h${}^{-1}$. The plausible range $[0.023,0.36]$ h${}^{-1}$ is guessed using the half-lives of N, M, E, and S proteins being about 30 and 1.9 h, respectively, in reticulocytes, and their relative molar ratio in a virion $\mathrm{E}$:$\mathrm{S}$:$\mathrm{M}=1$:20:300.

The composition of a single virion requires the binding of about 456 N protein molecules to each positive-sense genomic RNA, $n_N=456$. The total number of structural proteins $S,M,E$ can be estimated to be in between 1125 to 2230 [19,20,21] with the reference number we used about 2000, i.e., $n_{SP}=2000$.

The scale of SARS-CoV-2 replication depends on the target cell type (e.g., bronchial, lung cells, enterocytes). The available estimates suggest that the duration of the single replication cycle ranges from 7 to 24 h with the burst size in between 10 and 10,000 virions [2,22]. Hence, we set $K_{{\mathrm{V}}_{\mathrm{rel}}}=1000\in (10,10$,$000)$. Similar to modelling studies on replication of HIV-1 [10] and IAV [12], we estimate the following range for the rate of virion assembling $k_{assemb}$ = 0.01–10 h${}^{-1}$.

3.1.4. Assembly and Release of Virions

The assembly of virions requires that nucleocapsid and viral envelope glycoproteins coalesce into the same domain of the intracellular space [18]. The nucleocapsid core of the virion traffics to ERGIC and buds into ERGIC membranes covered with the structural proteins. Thus, a lipid envelope of the virion is created. The genome packaging is mediated by a packaging signal unique to genome length RNA.

As discussed above, the N proteins are key for incorporating viral RNA into viral progeny particles [38,39]. There, the N-terminal RNA-binding domain binds the RNA and the C-terminal domain via interaction with the M protein functions to anchor the ribonucleoprotein to the viral membrane [40].

The virions are assembled at the ER-Golgi compartment via encapsulating N-RNA complexes. Assembled new virions can exit the infected cell by exocytosis via the lysosomal trafficking pathway, budding, or cell death [17].

The rates of changes of the ribonucleocapsid and the assembled and released virions are described by the following equations:

(12)$\begin{array}{cc}\hfill \frac{d\left[\mathit{N-gRNA}\right]}{dt}& =k_{complex}{\theta }_{complex}\left[gRNA\right]-(k_{assemb}{\theta }_{assemb}+d_{N-gRNA})\left[\mathit{N-gRNA}\right]\hfill \end{array}$

(13)$\begin{array}{cc}\hfill \frac{d\left[V_{assembled}\right]}{dt}& =k_{assemb}{\theta }_{assemb}\left[\mathit{N-gRNA}\right]-(k_{release}+d_{assembled})\left[V_{assembled}\right]\hfill \end{array}$

(14)$\begin{array}{cc}\hfill \frac{d\left[V_{released}\right]}{dt}& =k_{release}\left[V_{assembled}\right]-d_V\left[V_{released}\right]\hfill \end{array}$

We use the estimate of $d_{gRNA}$ for $d_{N-gRNA}=0.2$ h${}^{-1}$, taking into account that the major component of ribonucleoprotein is ssRNA.

The direct measurements of the budding rate are not available yet. However, it has been shown for in vitro systems to be a very fast process with a characteristic time of about 2 s [41]. Hence, the following range is biologically plausible: $k_{release}\in [8,7200]$ h${}^{-1}$ [10].

As for the assembled virion death rate, we use the value $d_{assembled}=0.06$ h${}^{-1}$ estimated from [42], which is equal to $d_{endosome}$.

The overall list of model parameters with their reference values and permissible ranges is presented in Table 2. The solution of the model described by Equations (1)–(14) for $\left[V_{free}\right]\left(0\right)=10$ and the parameter values displayed in Table 2, is shown in Figure 2. It is consistent with data presented in [24,43].

In this study, we analyze the model behaviour at the initial condition $\left[V_{free}\right]\left(0\right)=10$. This corresponds to high MOI scenario of in vitro experiments in tissue cultures that are performed to estimate the burst size, i.e., the average number of virions produced by a single infected cell during the complete replication cycle. High MOI is typically used to ensure that every single cell gets infected and therefore only a single replication cycle occurs, resulting in a “one-step” growth dynamics of released progeny [2].

viruses-13-01735-t002_Table 2 Table 2

Estimates of the calibrated model parameters. The parameter range categories are labelled as follows: ranges based on the analysis of (†) literature, (‡) uncertainty quantification, ($†‡$) both.

Parameter	Description, Units	Value	Range, Relev. Refs.
$k_{bind}$	rate of virion binding to ACE2 receptor, h${}^{-1}$	12	$(3.6,12){}^{†}$ [44,45]
$d_V$	clearance rate of extracellular virions, h${}^{-1}$	$0.12$	$(0.06,3.5){}^{†}$ [42,46,47], tuned to $(0.06,0.2){}^{†‡}$
$k_{diss}$	dissociation rate constant of bound virions, h${}^{-1}$	$0.61$	$(0.32,1.08){}^{†}$ [44,45]
$k_{fuse}$	fusion rate constant, h${}^{-1}$	$0.5$	$(0.33,1){}^{†}$ [48]
$k_{uncoat}$	uncoating rate constant, h${}^{-1}$	$0.5$	$(0.33,1){}^{†}$ [48]
$d_{endosome}$	degradation rate of virions in endosomes, h${}^{-1}$	$0.06$	[12,42], $(0.0001,0.12){}^{‡}$
$k_{transl}$	translation rate, nt/mRNA h${}^{-1}$	45,360	[24,28], $(40$,$000,50$,$000){}^{‡}$
$1/f_{ORF1}$	length of ORF1 of the RNA genome coding $\left[NSP\right]$, nt	21,000	fixed [30]
$d_{NSP}$	degradation rate of proteins in the cell, h${}^{-1}$	$0.069$	$(0.023,0.69){}^{†}$ [28,37], tuned to $(0.023,0.1){}^{†‡}$
$k_{t{r}_{(-)}}$	transcription rate of negative sense genomic and subgenomic RNAs, copies/mRNA h${}^{-1}$	3	$(1,100){}^{†}$ [28], tuned to $(1,20){}^{†‡}$
$K_{NSP}$	threshold number of $\left[NSP\right]$ enhancing vRNA transcription, molecules	100	$(10,150){}^{‡}$
$d_{gRNA}$	degradation rate of positive sense RNAs in cell, h${}^{-1}$	$0.2$	$(0.069,0.69){}^{†}$ [28,31], tuned to $(0.069,0.4){}^{†‡}$
$d_{gRN{A}_{(-)}}$	degradation rate of negative sense RNAs in double-membrane vesicles, h${}^{-1}$	$0.1$	$(0.05,0.2){}^{‡}$
$k_{t{r}_{(+)}}$	replication rate of positive sense RNAs, copies/mRNA/h	1000	$(620,1380){}^{†}$ [32]
$k_{complex}$	rate of the nucleocapsid formation $[N-gRNA]$, h${}^{-1}$	$0.4$	$(0.02,0.4){}^{†}$ [21,33,49,50,51]
$K_N$	threshold number of N proteins at which nucleocapsid formation slows down, molecules	$5\times {10}^6$	$(3.5,6.5)\times {10}^{6†}$ [2,34,35]
$1/f_N$	length of RNA genome coding N protein, nt	1200	fixed [36]
$1/f_{SP}$	length of genome coding structural proteins $S,E,M$, nt	10,000	fixed [36]
$d_N$	degradation rate of N protein, h${}^{-1}$	$0.023$	$(0.023,0.069){}^{†}$ [37]
$d_{SP}$	mean degradation rate of the pool of $E,S,M$ proteins, h${}^{-1}$	$0.044$	$(0.023,0.36){}^{†}$ [37]
$n_{SP}$	total number of structural proteins $S,M,E$ per virion, molecules	2000	$(1125,2230){}^{†}$ [19,20,21]
$n_N$	number of N protein per virion, molecules	456	fixed [21]
$K_{{\mathrm{V}}_{\mathrm{rel}}}$	threshold number of virions at which the virion assembly process slows down, virions	1000	$(10,10$,$000){}^{†}$ [2,22]
$k_{assemb}$	rate of virion assembling, h${}^{-1}$	1	$(0.01,10){}^{†}$ [10,12]
$d_{\mathit{N-gRNA}}$	degradation rate of ribonucleoprotein, h${}^{-1}$	$0.2$	$(0.069,0.69){}^{†}$ [28,31]
$k_{release}$	rate of virion release via exocytosis, h${}^{-1}$	8	$(8,7200){}^{†}$ [10,41]
$d_{assembled}$	assembled virion degradation rate, h${}^{-1}$	$0.06$	[42], $(0.0001,0.12){}^{‡}$

The simulations cover the range of 24 h, which is informative for a single replication cycle setting experimental data on SARS-CoV-2 replication [2,22,23]. However, the interval can be stretched via amending the parameter values to adjust more specific data for particular cell types and infection conditions (e.g., Figure 3, right).

The model is validated by comparing its predictions against available experimental data as described in Methods section and shown in Figure 3.

3.2. Sensitivity Analysis

To evaluate the model response to parameter variations, a systematic sensitivity analysis of three model characteristics is performed. First, we assess the uncertainty of predicted progeny release kinetics caused by variations of the initial condition and model parameters in plausible ranges. Next, we identify the model parameters which have the most control over (1) the ratio of positive- to negative-sense viral RNAs and (2) the total number of virions released by 24 h post-infection.

3.2.1. Uncertainty Analysis of the Progeny Release Kinetics

As part of model validation, we assess the uncertainty of progeny release kinetics caused by two factors. First, we predicted the changes in progeny trajectories obtained with initial condition $\left[V_{free}\right]\left(0\right)$ varied by $\pm 50\%$, i.e., from 5 to 15 virions (Figure 4 (left)).

Next, we quantified the model uncertainty caused by the changes of model parameters (Figure 4 (right)). To this end, we randomly sampled a bunch of parameter combinations from the ranges of their permissible values as described in the Methods section. Model parameters presented in Table 2 can be either fixed ($f_{ORF1},f_N,f_{SP},n_N$) or allowed to be varied. For most of the later parameters, their permissible ranges were estimated from the literature. For parameters $k_{transl},d_{ensodome},K_{NSP},d_{gRN{A}_{(-)}},d_{assembled}$, we estimated their ranges so that the overall uncertainty of progeny release was matched with about 10 to 10,000 virions (i.e., 1000-fold range) in accordance with Table 1. However, we additionally needed to narrow down the ranges of some parameters which were determined from literature, i.e., $d_V,d_{NSP},k_{t{r}_{(-)}},d_{gRNA}$, so that the confidence region of progeny release at the end of replication cycle started from around ten virions but not from zero. As a result, this procedure allows us to consistently identify the plausible ranges for model parameters which ensure the robustness of the model behaviour in the range of its uncertainty that is matched to the level of uncertainty coming from current knowledge about the SARS-CoV-2 life cycle.

3.2.2. Parameters Controlling the Ratio of Positive- to Negative-Sense vRNAs

Given the good agreement of the calibrated model with available data on the kinetics of positive- to negative-sense vRNA ratio (Figure 3), we asked which model parameters determine this ratio. To this end, we analyze parameter sensitivity towards the following distance ${\mathrm{\Phi }}_{ratio}\left(\mathbf{p}\right)$ between the ratio predicted by the model $r_m\left(t_i\right)$ and ratios $r_{ex{p}_{\left(1\right)}}\left(t_i\right)$ and $r_{ex{p}_{\left(2\right)}}\left(t_i\right)$ from the two-replicate experiment [24]:

$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{ratio}\left(\mathbf{p}\right)=& \sum _i{\left({log}_{10}\left(r_m\left(t_i\right)\right)-{log}_{10}\left(r_{ex{p}_{\left(1\right)}}\left(t_i\right)\right)\right)}^2+\hfill \\ & \sum _i{\left({log}_{10}\left(r_m\left(t_i\right)\right)-{log}_{10}\left(r_{ex{p}_{\left(2\right)}}\left(t_i\right)\right)\right)}^2,r_m\left(t\right)=\frac{\left[gRN{A}_{(+)}\right]\left(t\right)+\left[gRNA\right]\left(t\right)}{\left[gRN{A}_{(-)}\right]\left(t\right)}.\hfill \end{array}$

The nonzero normalized sensitivity indices towards ${\mathrm{\Phi }}_{ratio}$ are shown in Figure 5. The following parameters have the largest effect on the discrepancy between the modelled and experimentally obtained ratio of positive- to negative-sense vRNA:

• threshold number of $\left[NSP\right]$ enhancing vRNA transcription,

• translation rate of non-structural proteins,

• rates of fusion and uncoating,

• replication rate of positive sense RNAs.

3.2.3. Predicting Novel Antiviral Targets That Control Progeny Production

The calibrated model can be used to predict the sensitivity of SARS-CoV-2 production by an infected cell to variations of the rates of underlying biochemical processes. To identify the prospective antiviral targets, we analyze parameter sensitivity towards virus progeny production. To this end, we use two functionals: (a) total number of released virions during $T=24$ h post-infection, ${\mathrm{\Phi }}_{progeny}\left(\mathbf{p}\right)={\int }_0^T{k}_{release}\left[V_{assembled}\right]dt\approx 280$ virions, and (b) the area under the curve type of metric for released virions during the same period, ${\mathrm{\Phi }}_{progeny}^{AUC}\left(\mathbf{p}\right)={\int }_0^T\left[V_{released}\right]dt\approx 1208$ virions × hours. The sensitivity values, normalized by parameter values, are presented in Figure 6 and Figure 7. The left part of the figures ranks the sensitivity indices for the model parameters that negatively impact the net SARS-CoV-2 production when their values are increased. The right bar plot ranks the sensitivity indices for the parameters which positively affect the net virus production with their increasing values. It follows that the three parameters having the largest sensitivity indices towards both functionals are the following ones:

• degradation rate of positive sense vRNAs in cytoplasm (negative effect),

• threshold number of $\left[NSP\right]$ enhancing vRNA transcription (negative effect),

• translation rate of non-structural proteins (positive effect).

4. Discussion

In this work, we have developed and calibrated a deterministic model of the SARS-CoV-2 life cycle at the level of the infected cell. The mathematical model describes the intracellular biochemistry underlying the replication of the virus. The model can be used as a part of multiscale models of the within-host SARS-CoV-2 infection, which are expected to provide a quantitative analytical tool to reveal complex pathogenetic mechanisms of COVID-19 following a systems approach [3].

The developed model predicts the dynamics of intracellular replication cycle of SARS-CoV-2 in target cells. The high initial condition of free infectious virions used in the paper corresponds to in vitro experimental studies of viral replication cycle with high MOI. This facilitates comparisons of model predictions of the dynamics of new viral progeny with those that would be estimated in future experimental setups. In such high MOI conditions, the entry of virions is very rapid [2]. In vivo, however, cells can be infected with a wide ranging number of free virions. Moreover, the expression of ACE2 on epithelial cells of various types and in different organs varies substantially [23,52], which would result in different kinetics of the virus entry. Therefore, the density of ACE2 on a cell membrane should be considered in the models of infection spreading within the human host organism. Furthermore, the lower the MOI, the more prevalent are stochastic effects at the early replication stages that result in the heterogeneity of the overall dynamics [53]. Thus, to predict the virus transmission in vivo, the multiscale hybrid models need to be developed that would incorporate the stochastic description of the intracellular SARS-CoV-2 life cycle.

Although many studies have been performed to characterize various aspects of the SARS-CoV-2 life cycle, detailed kinetics data similar to those available for HIV-1 [54] do not exist yet. Hence, the data “… to characterize the emergence of the viral replication intermediates and their impact on the cellular transcriptional response with high temporal resolution” [54] are urgently needed. In COVID-19, the virus infects multiple target cells expressing ACE2 including type I and II pneumocytes, alveolar macrophages, monocytes, endothelial cells, and airway epithelial cells, particularly of the mucous glands. To progress further with the development of relevant multiscale mathematical models, all these cell types need to be characterized with respect to the SARS-CoV-2 life cycle. For reviews on experimental models which are developed to study SARS-CoV-2 replication in cultures of various target cells, we refer to [55,56,57,58].

The deterministic quantitative description of the SARS-CoV-2 life cycle in the model establishes a mechanistic basis for the development of a stochastic description of the process kinetics using the Monte Carlo framework. This is required as many of the replication stages are characterized by low numbers of reactants and hence, the impact of random effects on the reaction kinetics is strong. The insight provided by stochastic models of influenza, hepatitis, and HIV infections [12,59] can be utilized in a similar way for SARS-CoV-2 modelling.

In our model, we do not consider spatial and transport aspects of the SARS-CoV-2 replication in the infected cell. An example of a stochastic approach to modelling the intracellular spatio-temporal dynamics of HIV-1 replication describing the microtubule transport of viral components is provided in [60]. In the case of SARS-CoV-2, adding the spatial dimension to the model would include the creation of replication compartments (e.g., double-membrane vesicles), intracellular trafficking (e.g., via the secretory pathway, and the molecular interactions with host proteins and innate immune responses [17]. For example, ORF3b, ORF6, and N proteins are known to interfere at multiple levels to inhibit the IFN-signalling pathway. The accumulation of viral factors leading to infected cell death via various mechanisms (apoptosis, necrosis, and pyroptosis) are not considered. These all require further analysis and strongly depend on the availability of appropriate quantitative data.

The developed mathematical model does not consider full details of subgenomic RNAs transcription and translation processes, and hence, the interaction of the encoded viral proteins with the intracellular defence mechanisms. The model needs to be linked to the cell physiology by considering the impact of virus infection on the codon usage, the cell metabolism, and the infected cell fate.

At present, there are no effective therapies against SARS-CoV-2 infection and computational modelling is used to assist with drug repurposing attempts [27]. Here, we have used our calibrated model to predict optimal targets for antiviral therapy. For this, all model parameters were ranked according to their contribution to the production of new virions by an infected cell using local sensitivity analysis. Parameters with the greatest negative and positive effects on virus progeny generation were identified. The three most influential parameters are (i) degradation rate of positive sense vRNAs in cytoplasm (negative effect), (ii) threshold number of non-structural proteins enhancing vRNA transcription (negative effect), and (iii) translation rate of non-structural proteins (positive effect). These parameters regulate the kinetics of various stages of the SARS-CoV-2 life cycle. Overall, systematic analyses of virus replication cycles with mathematical models of increasing level of detail should provide an efficient and rational way to define novel antiviral targets for therapies.

A quantitative mathematical description of the SARS-CoV-2 life cycle is a necessary step in gaining a predictive understanding of the regulatory mechanisms underlying the molecular biology of the virus and defining potential targets for inhibitors of the virus replication. Indeed, the sensitivity analysis of the model enables the prediction of potential targets for further experimental consideration, thus assisting the drug discovery or repurposing efforts. In fact, the derived model identifies such targets. The in silico model provides an analytical tool to examine the effects of certain virus mutations (e.g., the increase in affinity to ACE2 receptor) or combinations of multiple mutations on the viral progeny and drug resistance in advance of experiments. Finally, the model of SARS-CoV-2 replication being properly expanded in a question-driven manner could assist in gaining mechanistic insights into the robustness/fragility of intracellular defence mechanisms and the net outcomes of multi-modal therapeutic impacts which affect the specific stages of virus life cycle and the infected cell physiology.

Supplementary Materials

The following are available at https://www.mdpi.com/article/10.3390/v13091735/s1.

Author Contributions

Conceptualization, G.B., A.M. and D.G.; methodology, D.G., I.S. and G.B.; software, D.G.; validation, D.G. and I.S.; data curation, D.G., E.K. and I.S.; writing—review and editing, D.G., I.S., A.K., A.M. and G.B.; funding acquisition, G.B. All authors have read and agreed to the published version of the manuscript.

Funding

The reported study was funded by RFBR according to the research project numbers 20-04-60157 and 20-01-00352. AM is also supported by the Spanish Ministry of Science and Innovation grant no. PID2019-106323RB-I00(AEI/MINEICO/FEDER, UE), and “Unidad de Excelencia María de Maeztu”, funded by the AEI (CEX2018-000792-M). G.B. and D.G. were partly supported by Moscow Center for Fundamental and Applied Mathematics (agreement with the Ministry of Education and Science of the Russian Federation No. 075-15-2019-1624).

Data Availability Statement

Not applicable.

Acknowledgments

We thank the reviewers for their valuable suggestions on improving the clarity of the study.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

SARS-CoV-2	Severe acute respiratory syndrome coronavirus 2
ODE	Ordinary differential equations
MHV	Murine hepatitis virus
MOI	Multiplicity of infection

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Biochemical scheme of the SARS-CoV-2 replication cycle.

Figure 2. Reference model solution with parameters estimated in Table 2, [Vfree](0)=10. It predicts the kinetics of the viral replication intermediates. The upper row describes the cell entry state variables, i.e., (upper left) the free virions outside the cell membrane and the virions bound to ACE2, and (upper right) virions in endosomes together with single strand positive sense genomic RNA. The middle row specifies the kinetics of negative sense viral genome transcription and translation of non-structural proteins (middle left) followed by the positive sense genomic RNA transcription and the translation of N-protein (middle right). The bottom row displays the translation of the structural proteins and the ribonucleocapsid molecules formation (bottom left) resulting in creation of assembled virions in endosomes and the final release of virions (bottom right).

Figure 3. Validation of the calibrated model against available experimental data. (Left) ratio of positive- to negative-sense vRNA predicted by the model and measured in MHV-infected cells at MOI = 10 in [24]. (Right) SARS-CoV-2 progeny release kinetics in Vero E6 cells from [23,25] expressed in fold changes of the titer levels in cell culture from the start of their exponential growth. For model predictions, the moment of release onset was defined as the moment ts at which [Vreleased](ts)=1. The initial conditions were chosen to match the MOI used in experiments, i.e., so that the maximum number of uncoated [gRNA(+)]≈ MOI.

Figure 4. Uncertainty in the model output (progeny release kinetics). (Left) uncertainty associated with variation of the initial condition [Vfree](0). (Right) uncertainty associated with variation of model parameters in the ranges specified in Table 2.

Figure 5. Local normalized sensitivity indices having negative (left) and positive (right) effects on the objective function Φratio(p), i.e., on discrepancy between model and experiment.

Figure 6. Local normalized sensitivity indices having negative (left) and positive (right) effect on the total number of released virions Φprogeny(p).

Figure 7. Local normalized sensitivity indices having negative (left) and positive (right) effect on the area under the curve for released virions ΦprogenyAUC(p).