Next Article in Journal
Structural Analysis of the Menangle Virus P Protein Reveals a Soft Boundary between Ordered and Disordered Regions
Next Article in Special Issue
Personalized Virus Load Curves for Acute Viral Infections
Previous Article in Journal
Modified-Live Feline Calicivirus Vaccination Elicits Cellular Immunity against a Current Feline Calicivirus Field Strain in an Experimental Feline Challenge Study
Previous Article in Special Issue
Modeling Within-Host Dynamics of SARS-CoV-2 Infection: A Case Study in Ferrets
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences (INM RAS), 119333 Moscow, Russia
2
Moscow Center for Fundamental and Applied Mathematics at INM RAS, 119333 Moscow, Russia
3
World-Class Research Center “Digital Biodesign and Personalized Healthcare”, Sechenov First Moscow State Medical University, 119991 Moscow, Russia
4
Moscow Institute of Physics and Technology (National Research University), Dolgoprudny, 141701 Moscow Oblast, Russia
5
College of Engineering, Swansea University, Bay Campus, Fabian Way, Swansea SA1 8EN, UK
6
Department of Clinical Immunology and Allergology, Sechenov First Moscow State Medical University, 119991 Moscow, Russia
7
Infection Biology Laboratory, Universitat Pompeu Fabra, 08003 Barcelona, Spain
8
ICREA, Pg. Lluis Companys 23, 08010 Barcelona, Spain
9
Institute of Computer Science and Mathematical Modelling, Sechenov First Moscow State Medical University, 119991 Moscow, Russia
*
Authors to whom correspondence should be addressed.
Viruses 2021, 13(9), 1735; https://doi.org/10.3390/v13091735
Submission received: 21 May 2021 / Revised: 26 August 2021 / Accepted: 27 August 2021 / Published: 31 August 2021
(This article belongs to the Special Issue Mathematical Modeling of Viral Infection)

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.
Table 1. Time-dependent variables of the mathematical model characterizing the SARS-CoV-2 life cycle.
Table 1. Time-dependent variables of the mathematical model characterizing the SARS-CoV-2 life cycle.
VariableMeaningQuantitative Characteristics
[ V f r e e ] number of free virions outside the cell membrane10
[ V b o u n d ] number of virions bound to ACE2
and activated by TMPRSS2
1–10
[ V e n d o s o m e ] number of virions in endosomes1–10
[ g R N A ( + ) ] single strand positive sense genomic RNA1–5
[ N S P ] population of non-structural proteins
[ g R N A ( ) ] negative sense genomic and subgenomic RNAs10
[ g R N A ] positive sense genomic and subgenomic RNAs10,000
[ S P ] total number of structural proteins
S + M + E per virion
2000 ( 1125 , 2230 ) [19,20,21]
[ N ] N proteins per virion [ N ] 456 [21]; 1465 ( 730 , 2200 ) [19]
[ N-gRNA ] ribonucleocapsid molecules
[ V a s s e m b l e d ] assembled virions in endosomes
[ V r e l e a s e d ] virus burst size10–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 [ V r e l e a s e d ] ( t ) starting from the moment t s at which [ V r e l e a s e d ] ( t s ) = 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 [ V r e l e a s e d ] ( t ) trajectories were computed.

2.4. Sensitivity Analysis

Let u ( t , p ) be the model solution of the ODE initial value problem d u d t = f ( t , u , p ) , u ( 0 ) = 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 Φ ( u ( p ) ) , e.g., the cumulative number of released virions
Φ p r o g e n y = 0 T k r e l e a s e [ V a s s e m b l e d ] d t ,
the area under the curve for [ V r e l e a s e d ] ( t )
Φ p r o g e n y A U C ( p ) = 0 T [ V r e l e a s e d ] d t ,
or the distance between the model output g ( u ( p ) ) and the experimental data g e x p ( t i )
Φ = i   | | g ( u ( t i ) ) g e x p ( t i ) | | 2 .
The sensitivity index of model parameter p is defined as s p = d Φ ( p ) d p . It indicates the influence of parameter variation on the output value Φ . To compare sensitivities and rank parameters by their impact, the sensitivity indices are usually normalized by parameter values: s ^ p = p d Φ ( p ) d p . 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:
d [ V f r e e ] d t = k b i n d [ V f r e e ] d V [ V f r e e ] + k d i s s [ V b o u n d ]
d [ V b o u n d ] d t = k b i n d [ V f r e e ] ( k f u s e + k d i s s + d V ) [ V b o u n d ]
d [ V e n d o s o m e ] d t = k f u s e [ V b o u n d ] ( k u n c o a t + d e n d o s o m e ) [ V e n d o s o m e ]
d [ g R N A ( + ) ] d t = k u n c o a t [ V e n d o s o m e ] d g R N A [ g R N A ( + ) ] .
Here, [ V f r e e ] is the number of free virions outside the cell membrane, [ V b o u n d ] is the number of virions bound to ACE2 and activated by TMPRSS2, [ V e n d o s o m e ] is the number of virions in endosomes, and [ g R N A ( + ) ] 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 e n d o s o m e = 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 [ N S P ] , the set of negative sense genomic and subgenomic [ g R N A ( ) ] , and the set of positive sense genomic and subgenomic [ g R N A ] with the following differential equations:
d [ N S P ] d t = k t r a n s l f O R F 1 [ g R N A ( + ) ] d N S P [ N S P ]
d [ g R N A ( ) ] d t = k t r ( ) [ g R N A ( + ) ] θ R d R p d g R N A ( ) [ g R N A ( ) ]
d [ g R N A ] d t = k t r ( + ) [ g R N A ( ) ] θ R d R p ( k c o m p l e x θ c o m p l e x + d g R N A ) [ g R N A ]
where
θ R d R p = [ N S P ] [ N S P ] + K N S P , θ c o m p l e x = [ N ] [ N ] + 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 t r a n s l = 45 ,360 nt/mRNA h 1 has been estimated for coronaviruses in [24]. The length of the RNA genome coding for [ N S P ] proteins is about 21,000 nucleotides [30], hence f O R F 1 = 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 N S P = 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 N S P = 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 g R N A = 0.2 h 1 . In addition, we assume a smaller value for the degradation rate of negative sense vRNAs in double-membrane vesicles, d g R N A ( ) = d g R N A / 2 = 0.1 h 1 .
Quantitative analysis of RNA polymerase elongation provides the estimate of the transcription rate to be 46,080 ± 17,640 nt/h [32]. Thus, the basal rate of transcription is around 46,080 / 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 c o m p l e x . 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 c o m p l e x 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 × 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 μ M, and the concentration at which the formation slows down to be about 11 μ 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 × 10 6 molecules implies that the saturation takes place at the number of N proteins necessary to assemble around 11,000 virions ( K N n N · 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 [ S P ] . 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 [ N-g R N A ] by coating genomic RNAs. The number of N proteins per virion [ N ] 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 × 12 = 456 [21].
The translation rates of [ N ] and [ S P ] proteins are described by the following two equations:
d [ N ] d t = k t r a n s l f N [ g R N A ] k c o m p l e x n N θ c o m p l e x [ g R N A ] d N [ N ]
d [ S P ] d t = k t r a n s l f S P [ g R N A ] k a s s e m b n S P θ a s s e m b [ N-gRNA ] d S P [ S P ]
where
θ a s s e m b = [ S P ] [ S P ] + K V rel n S P
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 S P = 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 S P = 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 E : S : 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 S P = 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 V rel = 1000 ( 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 a s s e m b = 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:
d [ N-gRNA ] d t = k c o m p l e x θ c o m p l e x [ g R N A ] ( k a s s e m b θ a s s e m b + d N - g R N A ) [ N-gRNA ]
d [ V a s s e m b l e d ] d t = k a s s e m b θ a s s e m b [ N-gRNA ] ( k r e l e a s e + d a s s e m b l e d ) [ V a s s e m b l e d ]
d [ V r e l e a s e d ] d t = k r e l e a s e [ V a s s e m b l e d ] d V [ V r e l e a s e d ]
We use the estimate of d g R N A for d N - g R N A = 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 r e l e a s e [ 8 , 7200 ] h 1 [10].
As for the assembled virion death rate, we use the value d a s s e m b l e d = 0.06 h 1 estimated from [42], which is equal to d e n d o s o m e .
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 [ V f r e e ] ( 0 ) = 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 [ V f r e e ] ( 0 ) = 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].
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.
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.
ParameterDescription, UnitsValueRange, Relev. Refs.
k b i n d 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 d i s s dissociation rate constant of bound virions, h 1 0.61 ( 0.32 , 1.08 ) [44,45]
k f u s e fusion rate constant, h 1 0.5 ( 0.33 , 1 ) [48]
k u n c o a t uncoating rate constant, h 1 0.5 ( 0.33 , 1 ) [48]
d e n d o s o m e degradation rate of virions in endosomes, h 1 0.06 [12,42], ( 0.0001 , 0.12 )
k t r a n s l translation rate, nt/mRNA h 1 45,360[24,28],
( 40 , 000 , 50 , 000 )
1 / f O R F 1 length of ORF1 of the RNA genome coding [ N S P ] , nt21,000fixed [30]
d N S P 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 N S P threshold number of [ N S P ]
enhancing vRNA transcription, molecules
100 ( 10 , 150 )
d g R N A 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 g R N 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/h1000 ( 620 , 1380 ) [32]
k c o m p l e x rate of the nucleocapsid formation [ N - g R N A ] , 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 × 10 6 ( 3.5 , 6.5 ) × 10 6 [2,34,35]
1 / f N length of RNA genome coding N protein, nt1200fixed [36]
1 / f S P length of genome coding structural proteins S , E , M , nt10,000fixed [36]
d N degradation rate of N protein, h 1 0.023 ( 0.023 , 0.069 ) [37]
d S P mean degradation rate of the pool of E , S , M proteins, h 1 0.044 ( 0.023 , 0.36 ) [37]
n S P 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, molecules456fixed [21]
K V rel threshold number of virions
at which the virion assembly process slows down, virions
1000 ( 10 , 10 , 000 ) [2,22]
k a s s e m b rate of virion assembling, h 1 1 ( 0.01 , 10 ) [10,12]
d N-gRNA degradation rate of ribonucleoprotein, h 1 0.2 ( 0.069 , 0.69 ) [28,31]
k r e l e a s e rate of virion release via exocytosis, h 1 8 ( 8 , 7200 ) [10,41]
d a s s e m b l e d 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 [ V f r e e ] ( 0 ) varied by ± 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 O R F 1 , f N , f S P , n N ) or allowed to be varied. For most of the later parameters, their permissible ranges were estimated from the literature. For parameters k t r a n s l , d e n s o d o m e , K N S P , d g R N A ( ) , d a s s e m b l e d , 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 N S P , k t r ( ) , d g R N A , 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 Φ r a t i o ( p ) between the ratio predicted by the model r m ( t i ) and ratios r e x p ( 1 ) ( t i ) and r e x p ( 2 ) ( t i ) from the two-replicate experiment [24]:
Φ r a t i o ( p ) = i log 10 ( r m ( t i ) ) log 10 ( r e x p ( 1 ) ( t i ) ) 2 + i log 10 ( r m ( t i ) ) log 10 ( r e x p ( 2 ) ( t i ) ) 2 , r m ( t ) = [ g R N A ( + ) ] ( t ) + [ g R N A ] ( t ) [ g R N A ( ) ] ( t ) .
The nonzero normalized sensitivity indices towards Φ r a t i o 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 [ N S P ] 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, Φ p r o g e n y ( p ) = 0 T k r e l e a s e [ V a s s e m b l e d ] d t 280 virions, and (b) the area under the curve type of metric for released virions during the same period, Φ p r o g e n y A U C ( p ) = 0 T [ V r e l e a s e d ] d t 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 [ N S P ] 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

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-2Severe acute respiratory syndrome coronavirus 2
ODEOrdinary differential equations
MHVMurine hepatitis virus
MOIMultiplicity of infection

References

  1. Machhi, J.; Herskovitz, J.; Senan, A.M.; Dutta, D.; Nath, B.; Oleynikov, M.D.; Blomberg, W.R.; Meigs, D.D.; Hasan, M.; Patel, M.; et al. The Natural History, Pathobiology, and Clinical Manifestations of SARS-CoV-2 Infections. J. Neuroimmune Pharmacol. 2020, 15, 359–386. [Google Scholar] [CrossRef]
  2. Bar-On, Y.M.; Flamholz, A.; Phillips, R.; Milo, R. SARS-CoV-2 (COVID-19) by the numbers. eLife 2020, 9, e57309. [Google Scholar] [CrossRef]
  3. Sego, T.J.; Aponte-Serrano, J.O.; Ferrari Gianlupi, J.; Heaps, S.R.; Breithaupt, K.; Brusch, L.; Crawshaw, J.; Osborne, J.M.; Quardokus, E.M.; Plemper, R.K.; et al. A modular framework for multiscale, multicellular, spatiotemporal modeling of acute primary viral infection and immune response in epithelial tissues and its application to drug therapy timing and effectiveness. PLoS Comput. Biol. 2020, 16, e1008451. [Google Scholar] [CrossRef] [PubMed]
  4. Li, C.; Xu, J.; Liu, J.; Zhou, Y.; Li, C.; Xu, J.; Liu, J.; Zhou, Y. The within-host viral kinetics of SARS-CoV-2. Math. Biosci. Eng. 2020, 17, 2853–2861. [Google Scholar] [CrossRef] [PubMed]
  5. Wang, S.; Pan, Y.; Wang, Q.; Miao, H.; Brown, A.N.; Rong, L. Modeling the viral dynamics of SARS-CoV-2 infection. Math. Biosci. 2020, 328, 108438. [Google Scholar] [CrossRef]
  6. Hernandez-Vargas, E.A.; Velasco-Hernandez, J.X. In-host Mathematical Modelling of COVID-19 in Humans. Annu. Rev. Control 2020, 50, 448–456. [Google Scholar] [CrossRef]
  7. Liang, K. Mathematical model of infection kinetics and its analysis for COVID-19, SARS and MERS. Infect. Genet. Evol. 2020, 82, 104306. [Google Scholar] [CrossRef]
  8. Kim, K.S.; Ejima, K.; Iwanami, S.; Fujita, Y.; Ohashi, H.; Koizumi, Y.; Asai, Y.; Nakaoka, S.; Watashi, K.; Aihara, K.; et al. A quantitative model used to compare within-host SARS-CoV-2, MERS-CoV, and SARS-CoV dynamics provides insights into the pathogenesis and treatment of SARS-CoV-2. PLoS Biol. 2021, 19, e3001128. [Google Scholar] [CrossRef] [PubMed]
  9. Perelson, A.S.; Ke, R. Mechanistic Modeling of SARS-CoV-2 and Other Infectious Diseases and the Effects of Therapeutics. Clin. Pharmacol. Ther. 2021, 109, 829–840. [Google Scholar] [CrossRef] [PubMed]
  10. Shcherbatova, O.; Grebennikov, D.; Sazonov, I.; Meyerhans, A.; Bocharov, G. Modeling of the HIV-1 Life Cycle in Productively Infected Cells to Predict Novel Therapeutic Targets. Pathogens 2020, 9, 255. [Google Scholar] [CrossRef] [Green Version]
  11. Fatehi, F.; Bingham, R.J.; Dykeman, E.C.; Patel, N.; Stockley, P.G.; Twarock, R. An Intracellular Model of Hepatitis B Viral Infection: An In Silico Platform for Comparing Therapeutic Strategies. Viruses 2020, 13, 11. [Google Scholar] [CrossRef]
  12. Heldt, F.S.; Kupke, S.Y.; Dorl, S.; Reichl, U.; Frensing, T. Single-cell analysis and stochastic modelling unveil large cell-to-cell variability in influenza A virus infection. Nat. Commun. 2015, 6, 8938. [Google Scholar] [CrossRef]
  13. Aunins, T.R.; Marsh, K.A.; Subramanya, G.; Uprichard, S.L.; Perelson, A.S.; Chatterjee, A. Intracellular Hepatitis C Virus Modeling Predicts Infection Dynamics and Viral Protein Mechanisms. J. Virol. 2018, 92, e02098-17. [Google Scholar] [CrossRef] [Green Version]
  14. Teufel, A.I.; Liu, W.; Draghi, J.A.; Cameron, C.E.; Wilke, C.O. Modeling poliovirus replication dynamics from live time-lapse single-cell imaging data. Sci. Rep. 2021, 11, 9622. [Google Scholar] [CrossRef]
  15. Bocharov, G.; Casella, V.; Argilaguet, J.; Grebennikov, D.; Güerri-Fernandez, R.; Ludewig, B.; Meyerhans, A. Numbers Game and Immune Geography as Determinants of Coronavirus Pathogenicity. Front. Cell. Infect. Microbiol. 2020, 10, 559209. [Google Scholar] [CrossRef]
  16. Finkel, Y.; Gluck, A.; Nachshon, A.; Winkler, R.; Fisher, T.; Rozman, B.; Mizrahi, O.; Lubelsky, Y.; Zuckerman, B.; Slobodin, B.; et al. SARS-CoV-2 uses a multipronged strategy to impede host protein synthesis. Nature 2021, 594, 240–245. [Google Scholar] [CrossRef]
  17. V’kovski, P.; Kratzel, A.; Steiner, S.; Stalder, H.; Thiel, V. Coronavirus biology and replication: Implications for SARS-CoV-2. Nat. Rev. Microbiol. 2021, 19, 155–170. [Google Scholar] [CrossRef]
  18. Hartenian, E.; Nandakumar, D.; Lari, A.; Ly, M.; Tucker, J.M.; Glaunsinger, B.A. The molecular virology of coronaviruses. J. Biol. Chem. 2020, 295, 12910–12934. [Google Scholar] [CrossRef]
  19. Neuman, B.W.; Kiss, G.; Kunding, A.H.; Bhella, D.; Baksh, M.F.; Connelly, S.; Droese, B.; Klaus, J.P.; Makino, S.; Sawicki, S.G.; et al. A structural analysis of M protein in coronavirus assembly and morphology. J. Struct. Biol. 2011, 174, 11–22. [Google Scholar] [CrossRef]
  20. Yao, H.; Song, Y.; Chen, Y.; Wu, N.; Xu, J.; Sun, C.; Zhang, J.; Weng, T.; Zhang, Z.; Wu, Z.; et al. Molecular Architecture of the SARS-CoV-2 Virus. Cell 2020, 183, 730–738.e13. [Google Scholar] [CrossRef]
  21. Klein, S.; Cortese, M.; Winter, S.L.; Wachsmuth-Melm, M.; Neufeldt, C.J.; Cerikan, B.; Stanifer, M.L.; Boulant, S.; Bartenschlager, R.; Chlanda, P. SARS-CoV-2 structure and replication characterized by in situ cryo-electron tomography. Nat. Commun. 2020, 11, 5885. [Google Scholar] [CrossRef]
  22. Gordon, D.E.; Hiatt, J.; Bouhaddou, M.; Rezelj, V.V.; Ulferts, S.; Braberg, H.; Jureka, A.S.; Obernier, K.; Guo, J.Z.; Batra, J.; et al. Comparative host-coronavirus protein interaction networks reveal pan-viral disease mechanisms. Science 2020, 370, eabe9403. [Google Scholar] [CrossRef] [PubMed]
  23. Hou, Y.J.; Okuda, K.; Edwards, C.E.; Martinez, D.R.; Asakura, T.; Dinnon, K.H.; Kato, T.; Lee, R.E.; Yount, B.L.; Mascenik, T.M.; et al. SARS-CoV-2 Reverse Genetics Reveals a Variable Infection Gradient in the Respiratory Tract. Cell 2020, 182, 429–446.e14. [Google Scholar] [CrossRef]
  24. Irigoyen, N.; Firth, A.E.; Jones, J.D.; Chung, B.Y.W.; Siddell, S.G.; Brierley, I. High-Resolution Analysis of Coronavirus Gene Expression by RNA Sequencing and Ribosome Profiling. PLoS Pathog. 2016, 12, e1005473. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Ogando, N.S.; Dalebout, T.J.; Zevenhoven-Dobbe, J.C.; Limpens, R.W.; van der Meer, Y.; Caly, L.; Druce, J.; de Vries, J.J.C.; Kikkert, M.; Bárcena, M.; et al. SARS-coronavirus-2 replication in Vero E6 cells: Replication kinetics, rapid adaptation and cytopathology. J. Gen. Virol. 2020, 101, 925–940. [Google Scholar] [CrossRef] [PubMed]
  26. Marchuk, G.; Shutyaev, V.; Bocharov, G. Adjoint equations and analysis of complex systems: Application to virus infection modelling. J. Comput. Appl. Math. 2005, 184, 177–204. [Google Scholar] [CrossRef]
  27. Poduri, R.; Joshi, G.; Jagadeesh, G. Drugs targeting various stages of the SARS-CoV-2 life cycle: Exploring promising drugs for the treatment of Covid-19. Cell. Signal. 2020, 74, 109721. [Google Scholar] [CrossRef] [PubMed]
  28. Buccitelli, C.; Selbach, M. mRNAs, proteins and the emerging principles of gene expression control. Nat. Rev. Genet. 2020, 21, 630–644. [Google Scholar] [CrossRef]
  29. Sawicki, S.G.; Sawicki, D.L.; Siddell, S.G. A Contemporary View of Coronavirus Transcription. J. Virol. 2007, 81, 20–29. [Google Scholar] [CrossRef] [Green Version]
  30. Kim, D.; Lee, J.Y.; Yang, J.S.; Kim, J.W.; Kim, V.N.; Chang, H. The Architecture of SARS-CoV-2 Transcriptome. Cell 2020, 181, 914–921.e10. [Google Scholar] [CrossRef]
  31. Lehninger, A.L.; Nelson, D.L.; Cox, M.M. Lehninger Principles of Biochemistry, 5th ed.; W.H. Freeman: New York, NY, USA, 2008. [Google Scholar]
  32. Adelman, K.; La Porta, A.; Santangelo, T.J.; Lis, J.T.; Roberts, J.W.; Wang, M.D. Single molecule analysis of RNA polymerase elongation reveals uniform kinetic behavior. Proc. Natl. Acad. Sci. USA 2002, 99, 13538–13543. [Google Scholar] [CrossRef] [Green Version]
  33. Zinzula, L.; Basquin, J.; Bohn, S.; Beck, F.; Klumpe, S.; Pfeifer, G.; Nagy, I.; Bracher, A.; Hartl, F.U.; Baumeister, W. High-resolution structure and biophysical characterization of the nucleocapsid phosphoprotein dimerization domain from the Covid-19 severe acute respiratory syndrome coronavirus 2. Biochem. Biophys. Res. Commun. 2021, 538, 54–62. [Google Scholar] [CrossRef] [PubMed]
  34. Jack, A.; Ferro, L.S.; Trnka, M.J.; Wehri, E.; Nadgir, A.; Nguyenla, X.; Costa, K.; Stanley, S.; Schaletzky, J.; Yildiz, A. SARS-CoV-2 nucleocapsid protein forms condensates with viral genomic RNA. bioRxiv 2021. [Google Scholar] [CrossRef]
  35. Cubuk, J.; Alston, J.J.; Incicco, J.J.; Singh, S.; Stuchell-Brereton, M.D.; Ward, M.D.; Zimmerman, M.I.; Vithani, N.; Griffith, D.; Wagoner, J.A.; et al. The SARS-CoV-2 nucleocapsid protein is dynamic, disordered, and phase separates with RNA. Nat. Commun. 2021, 12, 1936. [Google Scholar] [CrossRef] [PubMed]
  36. Viehweger, A.; Krautwurst, S.; Lamkiewicz, K.; Madhugiri, R.; Ziebuhr, J.; Hölzer, M.; Marz, M. Direct RNA nanopore sequencing of full-length coronavirus genomes provides novel insights into structural variants and enables modification analysis. Genome Res. 2019, 29, 1545–1554. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Gasteiger, E.; Gattiker, A.; Hoogland, C.; Ivanyi, I.; Appel, R.D.; Bairoch, A. ExPASy: The proteomics server for in-depth protein knowledge and analysis. Nucleic Acids Res. 2003, 31, 3784–3788. [Google Scholar] [CrossRef] [Green Version]
  38. Khan, M.T.; Irfan, M.; Ahsan, H.; Ahmed, A.; Kaushik, A.C.; Khan, A.S.; Chinnasamy, S.; Ali, A.; Wei, D.Q. Structures of SARS-CoV-2 RNA-Binding Proteins and Therapeutic Targets. Intervirology 2021, 64, 55–68. [Google Scholar] [CrossRef]
  39. Khan, M.T.; Zeb, M.T.; Ahsan, H.; Ahmed, A.; Ali, A.; Akhtar, K.; Malik, S.I.; Cui, Z.; Ali, S.; Khan, A.S.; et al. SARS-CoV-2 nucleocapsid and Nsp3 binding: An in silico study. Arch. Microbiol. 2021, 203, 59–66. [Google Scholar] [CrossRef]
  40. Surjit, M.; Lal, S.K. The Nucleocapsid Protein of the SARS Coronavirus: Structure, Function and Therapeutic Potential. In Molecular Biology of the SARS-Coronavirus; Lal, S.K., Ed.; Springer: Berlin/Heidelberg, Germany, 2010; pp. 129–151. [Google Scholar] [CrossRef]
  41. Mooney, J.; Thakur, S.; Kahng, P.; Trapani, J.G.; Poccia, D. Quantification of exocytosis kinetics by DIC image analysis of cortical lawns. J. Chem. Biol. 2014, 7, 43–55. [Google Scholar] [CrossRef] [Green Version]
  42. Baggen, J.; Persoons, L.; Vanstreels, E.; Jansen, S.; Van Looveren, D.; Boeckx, B.; Geudens, V.; De Man, J.; Jochmans, D.; Wauters, J.; et al. Genome-wide CRISPR screening identifies TMEM106B as a proviral host factor for SARS-CoV-2. Nat. Genet. 2021, 53, 435–444. [Google Scholar] [CrossRef]
  43. Lokugamage, K.G.; Hage, A.; de Vries, M.; Valero-Jimenez, A.M.; Schindewolf, C.; Dittmann, M.; Rajsbaum, R.; Menachery, V.D. Type I Interferon Susceptibility Distinguishes SARS-CoV-2 from SARS-CoV. J. Virol. 2020, 94, e01410-20. [Google Scholar] [CrossRef]
  44. Ozono, S.; Zhang, Y.; Ode, H.; Sano, K.; Tan, T.S.; Imai, K.; Miyoshi, K.; Kishigami, S.; Ueno, T.; Iwatani, Y.; et al. SARS-CoV-2 D614G spike mutation increases entry efficiency with enhanced ACE2-binding affinity. Nat. Commun. 2021, 12, 848. [Google Scholar] [CrossRef] [PubMed]
  45. Walls, A.C.; Park, Y.J.; Tortorici, M.A.; Wall, A.; McGuire, A.T.; Veesler, D. Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein. Cell 2020, 181, 281–292.e6. [Google Scholar] [CrossRef] [PubMed]
  46. Bocharov, G.; Romanyukha, A. Mathematical Model of Antiviral Immune Response III. Influenza A Virus Infection. J. Theor. Biol. 1994, 167, 323–360. [Google Scholar] [CrossRef] [PubMed]
  47. Baccam, P.; Beauchemin, C.; Macken, C.A.; Hayden, F.G.; Perelson, A.S. Kinetics of Influenza A Virus Infection in Humans. J. Virol. 2006, 80, 7590–7599. [Google Scholar] [CrossRef] [Green Version]
  48. Zhu, Y.; Yu, D.; Yan, H.; Chong, H.; He, Y. Design of Potent Membrane Fusion Inhibitors against SARS-CoV-2, an Emerging Coronavirus with High Fusogenic Activity. J. Virol. 2020, 94, e00635-20. [Google Scholar] [CrossRef]
  49. Chen, I.J.; Yuann, J.M.P.; Chang, Y.M.; Lin, S.Y.; Zhao, J.; Perlman, S.; Shen, Y.Y.; Huang, T.H.; Hou, M.H. Crystal structure-based exploration of the important role of Arg106 in the RNA-binding domain of human coronavirus OC43 nucleocapsid protein. Biochim. Biophys. Acta (BBA)-Proteins Proteom. 2013, 1834, 1054–1062. [Google Scholar] [CrossRef]
  50. Spencer, K.A.; Hiscox, J.A. Characterisation of the RNA binding properties of the coronavirus infectious bronchitis virus nucleocapsid protein amino-terminal region. FEBS Lett. 2006, 580, 5993–5998. [Google Scholar] [CrossRef] [Green Version]
  51. Spencer, K.A.; Dee, M.; Britton, P.; Hiscox, J.A. Role of phosphorylation clusters in the biology of the coronavirus infectious bronchitis virus nucleocapsid protein. Virology 2008, 370, 373–381. [Google Scholar] [CrossRef] [Green Version]
  52. Ziegler, C.G.; Allon, S.J.; Nyquist, S.K.; Mbano, I.M.; Miao, V.N.; Tzouanas, C.N.; Cao, Y.; Yousif, A.S.; Bals, J.; Hauser, B.M.; et al. SARS-CoV-2 Receptor ACE2 Is an Interferon-Stimulated Gene in Human Airway Epithelial Cells and Is Detected in Specific Cell Subsets across Tissues. Cell 2020, 181, 1016–1035.e19. [Google Scholar] [CrossRef]
  53. Sazonov, I.; Grebennikov, D.; Kelbert, M.; Meyerhans, A.; Bocharov, G. Viral Infection Dynamics Model Based on a Markov Process with Time Delay between Cell Infection and Progeny Production. Mathematics 2020, 8, 1207. [Google Scholar] [CrossRef]
  54. Mohammadi, P.; Desfarges, S.; Bartha, I.; Joos, B.; Zangger, N.; Muñoz, M.; Günthard, H.F.; Beerenwinkel, N.; Telenti, A.; Ciuffi, A. 24 Hours in the Life of HIV-1 in a T Cell Line. PLoS Pathog. 2013, 9, e1003161. [Google Scholar] [CrossRef] [PubMed]
  55. De Dios-Figueroa, G.T.; Aguilera-Marquez, J.d.R.; Camacho-Villegas, T.A.; Lugo-Fabres, P.H. 3D Cell Culture Models in COVID-19 Times: A Review of 3D Technologies to Understand and Accelerate Therapeutic Drug Discovery. Biomedicines 2021, 9, 602. [Google Scholar] [CrossRef]
  56. Heinen, N.; Klöhn, M.; Steinmann, E.; Pfaender, S. In Vitro Lung Models and Their Application to Study SARS-CoV-2 Pathogenesis and Disease. Viruses 2021, 13, 792. [Google Scholar] [CrossRef]
  57. Synowiec, A.; Szczepański, A.; Barreto-Duran, E.; Lie, L.K.; Pyrc, K. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2): A Systemic Infection. Clin. Microbiol. Rev. 2021, 34, e00133-20. [Google Scholar] [CrossRef]
  58. Simoneau, C.R.; Ott, M. Modeling Multi-organ Infection by SARS-CoV-2 Using Stem Cell Technology. Cell Stem Cell 2020, 27, 859–868. [Google Scholar] [CrossRef] [PubMed]
  59. Andreu-Moreno, I.; Bou, J.V.; Sanjuán, R. Cooperative nature of viral replication. Sci. Adv. 2020, 6, eabd4942. [Google Scholar] [CrossRef]
  60. Cheng, Z.; Hoffmann, A. A stochastic spatio-temporal (SST) model to study cell-to-cell variability in HIV-1 infection. J. Theor. Biol. 2016, 395, 87–96. [Google Scholar] [CrossRef]
Figure 1. Biochemical scheme of the SARS-CoV-2 replication cycle.
Figure 1. Biochemical scheme of the SARS-CoV-2 replication cycle.
Viruses 13 01735 g001
Figure 2. Reference model solution with parameters estimated in Table 2, [ V f r e e ] ( 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 2. Reference model solution with parameters estimated in Table 2, [ V f r e e ] ( 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).
Viruses 13 01735 g002
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 t s at which [ V r e l e a s e d ] ( t s ) = 1 . The initial conditions were chosen to match the MOI used in experiments, i.e., so that the maximum number of uncoated [ g R N A ( + ) ] MOI.
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 t s at which [ V r e l e a s e d ] ( t s ) = 1 . The initial conditions were chosen to match the MOI used in experiments, i.e., so that the maximum number of uncoated [ g R N A ( + ) ] MOI.
Viruses 13 01735 g003
Figure 4. Uncertainty in the model output (progeny release kinetics). (Left) uncertainty associated with variation of the initial condition [ V f r e e ] ( 0 ) . (Right) uncertainty associated with variation of model parameters in the ranges specified in Table 2.
Figure 4. Uncertainty in the model output (progeny release kinetics). (Left) uncertainty associated with variation of the initial condition [ V f r e e ] ( 0 ) . (Right) uncertainty associated with variation of model parameters in the ranges specified in Table 2.
Viruses 13 01735 g004
Figure 5. Local normalized sensitivity indices having negative (left) and positive (right) effects on the objective function Φ r a t i o ( p ) , i.e., on discrepancy between model and experiment.
Figure 5. Local normalized sensitivity indices having negative (left) and positive (right) effects on the objective function Φ r a t i o ( p ) , i.e., on discrepancy between model and experiment.
Viruses 13 01735 g005
Figure 6. Local normalized sensitivity indices having negative (left) and positive (right) effect on the total number of released virions Φ p r o g e n y ( p ) .
Figure 6. Local normalized sensitivity indices having negative (left) and positive (right) effect on the total number of released virions Φ p r o g e n y ( p ) .
Viruses 13 01735 g006
Figure 7. Local normalized sensitivity indices having negative (left) and positive (right) effect on the area under the curve for released virions Φ p r o g e n y A U C ( p ) .
Figure 7. Local normalized sensitivity indices having negative (left) and positive (right) effect on the area under the curve for released virions Φ p r o g e n y A U C ( p ) .
Viruses 13 01735 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Grebennikov, D.; Kholodareva, E.; Sazonov, I.; Karsonova, A.; Meyerhans, A.; Bocharov, G. Intracellular Life Cycle Kinetics of SARS-CoV-2 Predicted Using Mathematical Modelling. Viruses 2021, 13, 1735. https://doi.org/10.3390/v13091735

AMA Style

Grebennikov D, Kholodareva E, Sazonov I, Karsonova A, Meyerhans A, Bocharov G. Intracellular Life Cycle Kinetics of SARS-CoV-2 Predicted Using Mathematical Modelling. Viruses. 2021; 13(9):1735. https://doi.org/10.3390/v13091735

Chicago/Turabian Style

Grebennikov, Dmitry, Ekaterina Kholodareva, Igor Sazonov, Antonina Karsonova, Andreas Meyerhans, and Gennady Bocharov. 2021. "Intracellular Life Cycle Kinetics of SARS-CoV-2 Predicted Using Mathematical Modelling" Viruses 13, no. 9: 1735. https://doi.org/10.3390/v13091735

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop