Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) virus is an intracellular infectious agent whose replication cycle depends on the host’s cell structures and functions. In particular, it uses the translational apparatus of different types of infected cells to express its proteins (Nakagawa et al., 2016). SARS-CoV-2 causes the Coronavirus Disease 2019 (COVID19) that has infected over 241, 000, 000 persons and killed over 4,900,000 worldwide since the end of 2019. No specific and effective therapeutic drugs to defeat SARS-CoV-2 infection have been developed yet (Díaz, 2020a), although several effective vaccines have been developed and are being applied. Therefore, it is necessary to continue the search for effective therapeutic agents that complement the use of vaccines, to minimize the risk of death from COVID19.

SARS-CoV-2 virion is formed by four proteins: spike (S), envelope (E), membrane (M) and nucleocapside (N) that enclose the virus genome (McBride and Fielding, 2012), which consists of a positive-sense nonsegmented single stranded mRNA [(+)ssRNA] of 30 kb. The open reading frames 1a (orf1a) and 1b (orf1b) are located near the untranslated 5′ region (5′UTR) of the positive single stranded RNA [(+)ssRNA] and they encode for the polyproteins pp1a and pp1ab.

When SARS-CoV-2 virion infects the organism, S protein binds with high affinity to the surface receptor Angiotensin-Converting Enzyme 2 (ACE2), highly abundant in the lung alveolar type II cells and other cell types (Hamming et al., 2004; Gordon et al., 2020), and forms a molecular complex that begins the process of fusion of the virion envelope with the host cell membrane. Finally, viral (+)ssRNA is released into the host cytoplasm (Letko et al., 2020) (Figure 1).

The process of translation occurs in the cytoplasm and produces a set of viral polyproteins whose maturation results in 11 nonstructural proteins (Nsp) from the orf1a segment (Nsp1 to Nsp11) and 5 nonstructural proteins from the orf1b segment (Nsp12 to Nsp16). Nsp proteins form the replication-transcription complex (RTC) in a double-membrane vesicle where a set of nested subgenomic minus-strands of RNA ((-)sgRNA) are synthesized in a process of discontinuous transcription. These (-) sgRNAs serve as the templates for the production of subgenomic mRNAs from which the structural proteins E, M, N and S, together with the accessory proteins orf 3a, orf6, orf7a, orf7b, orf8, orf9b, orf9c and orf10 are synthesized (Sevajol et al., 2014; Díaz, 2020a; Dongwan et al., 2020). SARS-CoV-2 uses the host translational machinery to redirect it to viral protein synthesis and replication, while host mRNA translation is inhibited (Nakagawa et al., 2016). SARS-CoV and SARS-CoV-2 genomes have ∼79% of homology (Forster et al., 2020), and most of the set of structural and nonstructural proteins are practically the same. However, the virus species differ in the accessory proteins orf8, orf8a, orf 8b, orf9c and orf10 (Bartlam et al., 2005; Gordon et al., 2020).

Viral proteins are inserted into the host molecular machinery to modify and redirect a great number of host cell functions towards the production of additional viral particles (Masters, 2006; Wu et al., 2020). Experimental analysis of the interaction of viral and host proteins, or interactome, by Gordon and collaborators (Gordon et al., 2020) has been a fundamental contribution to understand the form in which SARS-CoV-2 virus takes control of the host molecular network to produce new virions and propagate the infection (Hekman et al., 2020). Gordon and collaborators (2020) cloned, tagged and expressed 26 viral proteins in human cells using affinity-purification mass spectrometry to identify the human proteins physically associated with each other. They found around 332 SARS-CoV-2-human protein-protein interactions that constitute the virus interactome. From these results, the construction of the network representation of the interactome is possible. In the particular case of the Gordon interactome, the analysis of its undirected network model indicates a modular free-scale hierarchical type of structure in which proteins orf8, N and Nsp7 are the main hubs (Díaz, 2020a).

Despite all the molecular knowledge rapidly accumulated on SARS-CoV-2 itself and the infection or inflammation processes in which it is involved in human cells, we still lack dynamical models to understand which are the key nodes underlying severe inflammation processes in some cases. Mathematical and computational dynamic network models are useful tools to integrate data, find experimental holes and also propose or understand the underlying dynamics mechanisms involved in the complex viral/human networks involved during SARS-CoV-2 infection (Breitling, 2010; Díaz, 2020a; Enciso et al., 2020; Weinstein et al., 2020). The number of nodes, the number of connections of each node to its neighbours, and the distribution of these connections in the network determine its complexity, structure and dynamical properties (Kumar et al., 2020).

In this network, the viral main protease Nsp5 is the main protease (Mpro) involved in the assembly of the polyprotein cluster that translates viral RNA, and in the post-translational modification of viral proteins through its ADP ribose phosphatase activity (Nakagawa et al., 2016). SARS-CoV-2 Nsp5 has 96% sequence homology with its SARS-CoV counterpart, and its level of phylogenetic conservation by sequence and structure is high having high homogeneity in its volume and electrostatic profile throughout the Coronaviridae family. Interest in elucidating the functional structure of this protease has resulted in new reports of 3D structures showing that Nsp5 is a complex formed by two equal protomers (A and B) each consisting of 306 residues structured in three domains: I Chymotrypsin-like domain (residues 8–101), Picornavirus 3C protease-like domain II (residues 102–184), and Globular Cluster domain III (residues 201–303). Domains I and II have an antiparallel β-barrel structure and associate with their counterparts to form a dimer, which loop residues 185–200 connect to domain III made up of five α-helixes. The substrate binding site is found in the antiparallel β-barrels located between domains I and II, forming a catalytic pocket composed of four subsites (S1, S2, S3 and S4) that exhibit charged residues that include the nucleophilic sulphur atom of Cys145, which forms a dyad with the imidazole ring of His 41. Therefore, the blocking one of these residues with a small molecule could inhibit the enzymatic activity of Nsp5 (Mukherjee et al., 2011; Pillaiyar et al., 2016; Konwar and Sarma, 2021).

There are several crystallographic structures of Mpro with peptidomimetic inhibitors that show interactions with the enzyme at S3-S4 subsites in SARS-CoV-2, SARS-CoV and MERS. This type of inhibitor incorporates an imidazole ring of cyclic ketoamides with five-components that binds covalently to Cys145 or its equivalent in the different isoforms of the enzyme. Although these types of molecules are not optimal for drug development, due to proteolytic degradation, they have made possible to describe the dynamics of the enzymatic mechanism of this active site for the search and development of molecules with better pharmacodynamic parameters. An example of viable candidates for clinical implementation is boceprevir and telaprevir, whose mechanism of action consists of blocking the catalytic dyad through the covalent binding to the sulphur atom of Cys145 (Pan et al., 2012).

In the SARS-CoV-2 molecular network, Nsp5 has a special role in lung inflammation process through its link with Histone Deacetylase 2 (HDAC2) (Díaz, 2020a; Gordon et al., 2020; Hekman et al., 2020) (Figures 1, 2); suggesting that Nsp5 represses HDAC2 inhibitory action on NF κB (Figure 1), and thus enhances the transcription of the pro-inflammatory genes targeted by the Nuclear Factor κB (NF κB) (Figures 1, 2) (Wagner et al., 2015). Such role distinguishes SARS-CoV and SARS-CoV-2. The SARS-CoV virus uses N protein to promote the sustained transcription of the Cyclooxygenase 2 (Cox2) enzyme by binding directly to the NF-κB transcriptional regulatory elements, and to the CCAAT/enhancer binding proteins of the gene Cox2 (Yan et al., 2006). Hence, SARS-CoV and SARS-CoV-2 produce severe acute respiratory syndrome through different molecular mechanisms.

The change from a hub (N protein in SARS-CoV) to a poor connected protein (Nsp5 in SARS-CoV-2) as the underlying cause of severe acute respiratory syndrome could be an adaptation that made SARS-CoV-2 more pathogenic (8,098 people worldwide became sick with SARS during the 2003 outbreak against 241, 000, 000 during the 2020–2021 SARS-Cov-2 outbreak). Nsp5 is also the main protease that cleavages the SARS-CoV-2 polyprotein at 11 conserved sites, and its malfunction or deletion could stop the viral replication cycle and lung inflammation (Cohen et al., 2000; Xun et al., 2020). Thus, Nsp5 is a molecular switch that shows a low rate of mutation (Vilar and Isom, 2020), and its malfunctioning can decrease the strength of inflammation and the persistence of ACE2 in the lung epithelial cells giving rise to a less intense inflammatory response (Cohen et al., 2000; Xun et al., 2020).

Activation of the NF κB family of proteins (p65/RelA, p105/p50, p100/p52, RelB, and c-Rel) is mediated by the phosphorylation of the inhibitor of NF κB proteins (IκBs) by IκB kinases (IKKs) and the posterior ubiquitinylation and degradation of the phosphorylated IκBs (Figure 1). The released NF κB proteins enter the nucleus to activate specific genes (Wagner et al., 2015). In the nucleus, HDAC2 blocks the transcriptional activity of NF κB, inhibiting the production and release of the cytokine Interleukin 6 (IL-6) and IκB, and interfering with the functioning of the positive feedback circuit between IL-6 and NF κB (Rahman, and MacNee, 1998) (Figure 3A). As we mentioned before, the sustained inactivation of HDAC2 by Nsp5 produces a sustained increase in the transcriptional activity of NF κB (Figures 1, 2) that leads to a high production of IL-6 in the lung epithelium cells. Release of this excess of IL-6 to the vascular system contributes to the “cytokines storm” observed in critical patients (Leisman et al., 2020; Magro, 2020), although the average serum concentration of IL-6 in critical COVID19 patients are lower with respect to other respiratory syndromes (Leisman et al., 2020). However, the uninterrupted production of IL-6 is closely associated to a persistent presence of the ACE2 receptor in the alveolar type 2 cells due to enhanced ACE2 gene expression mediated by the JAK-STAT and PI3K/Akt pathways (Zegeye et al., 2018; Hennighausen and Lee, 2020; Mokuda et al., 2020) (Figure 2).

NF κB also promotes transcription of the Cox2 gene leading to the expression and increased enzymatic activity of Cox2 enzyme (Figures 2, 3B). However, NF κB, Cox2 and IL-6 also form a positive feedback circuit in which NF κB promotes Cox2 enzymatic activity that catalyses the conversion of Arachidonic acid (AA) into Prostaglandin G2 (PGG2), which is in turn modified by the peroxidase moiety of the Cox2 enzyme to produce prostaglandin H2 (PGH2) that is converted to prostaglandin E2 (PGE2) (Alexanian and Sorokin, 2017). PGE2 is then released to the vascular system, binds to its membrane EP2 and EP4 prostanoid receptors and promotes the cleavage of the IκB-NF κB complex increasing the amount of free NF κB, which increases the production of IL-6 (Figure 3B) (Cho et al., 2014; Bouffi et al., 2010). IL-6, in turn, increases NF κB activity (Wang et al., 2003) (Figure 2). These regulatory circuits has three therapeutic targets: IL-6 production (Tocilizumab) (Magro 2020), Cox2 (Nimesulide) (Alexanian and Sorokin, 2017) and NF κB (Dexamethasone) (Newton et al., 1998; Aghai et al., 2006). However, none of these have shown to be really effective for preventing or reverting acute inflammatory cases after SARS-CoV-2 infection.

In this work, we propose a dynamical mathematical model that suggests that the three feedback loops mentioned above constitute a minimal circuit for the inflammatory process in the epithelial lung cells during SARS-CoV-2 infection (Figure 2), and we use an ordinary differential equations (ODEs) continuous model to explore the effect of Nsp5 in the qualitative dynamics of this circuit in which this viral protein acts as an enhancing perturbation in the phase space of this dynamical system (Díaz, 2020b) (Figure 2). For the model of the circuit, we chose as main molecular components the proteins NF κB, IκB, IL-6, Cox2, and PGE2. In this model, Nsp5 prevents the movement of HDAC2 into the nucleus allowing an increase in NF κB transcriptional activity. The input that turns on the circuit in lung cells is a signal from monocytes, lymphocytes, and other immune cells that activates IL-6 production. We propose that this model represents the basic set of interactions that settles on the inflammation reaction to the invasion of the lung cells by SARS-CoV-2, and that it is a first approximation to understand the main dynamical features of this response.

Our hypothesis is that Nsp5 leads to the sustained overproduction of IL-6, Cox2, PGE2, and NF κB, in the epithelial lung cells, which is a necessary but not a sufficient condition for a cytokine storm. In the particular case of epithelial lung cells, we found that the computational and mathematical analysis of the model supports our hypothesis, and that Nsp5 effectively enhances NF κB, IL-6, Cox2 and PGE2 production during the time that remains bounded to HDAC2. Furthermore, the model predicts the existence of an over stimulated immune state (OSIS) in which the amount of IL-6 produced by the infected lung cells is high, and such state becomes autonomous and robust, independent of any stimulus from external sources like monocytes, and lymphocytes. Hence this OSIS cannot be completely shut down by classical anti-inflammatory drugs like Nimesulide and Dexamethasone, and can be possibly beaten only with specific inhibitors of Nsp5 (Xun et al., 2020).

Immune response to infections involves complex spatio-temporal processes involving multiple molecular components, non-linear interactions and cell-cell interactions. In consequence, the immune response to infections is a multidimensional, multiscale and multilevel dynamical process in which innate and acquired immunity mechanisms are involved to protect the organism against pathogens (Eftimie et al., 2016). In this scenario, the number of variables that become both sufficient and necessary for a detailed enough description of the spatio-temporal dynamics of the immune system is tremendously high. Some computational integrative efforts have been made in this direction (Danos et al., 2007), but imply very large systems that cannot be modeled with continuous quantitative models that are, in turn, necessary if we want to understand which are key nodes driving the system towards the OSIS. The lack of quantitative experimental data and knowledge of the values of most of the model parameters are additional challenges to achieve such dynamic models of the immune system in response to specific infectious diseases. Thus, most models of the immune system have been qualitative (Eftimie et al., 2016; Martinez-Sanchez et al., 2018).

A possible approach implies the use of ordinary differential equations (ODEs) qualitative modeling that is oriented to determine all the possible trajectories of a dynamical system along its n-dimensional phase space (with n > 0), when the system is subject to a set of initial conditions and parameters values. High dimensional nonlinear dynamical systems can exhibit a variety of coexisting dynamical behaviors that determine the structure of their phase space (Strogatz, 2015). However, the structure of the phase space of biological nonlinear dynamical systems is highly dependent on the set of parameters values. Variations in the value of one or more parameters can drive a drastic change in phase space generating a different qualitative dynamics of the system. This process is known as bifurcation, and it is important to find the possible bifurcations in the system and to identify which parameter or parameters are the responsible of these changes (Strogatz, 2015) (Supplementary Material S1).

In order to avoid the problem due to the high dimensionality and complexity of the immune response to SARS-CoV-2 infection, and to use the minimum number of parameters in the model that we are proposing in this work, we assume that the interaction between nodes IL-6, NF κB, IκB, Cox2 and PGE2 settles on the main dynamical features of the inflammation process in lung epithelial cells due to the presence of the Nsp5 viral protein. In the proposed circuit shown in Figure 2, IL-6, NF κB, Cox2 and PGE2 define three feedback loops (Figures 3A,B): a positive feedback loop between NF κB and IL-6 (Figure 3A); a positive feedback loop between NF κB and IκB (Figure 3A); and a negative feedback loop between PGE2 and IκB (Figure 3B).

These feedback loops are distributed between three compartments:

1) IL-6 is secreted into the extracellular medium where it acts in paracrine and exocrine form. IL-6 binds to its receptor IL6R at the cell surface. Once the complex IL6-IL6R is formed, it activates IKK in the cytoplasm (Wang et al., 2003). PGE2 is also secreted into this compartment where it binds to EP2 and EP4 receptors at the cell surface, activating Protein Kinase B (Akt) in the cytoplasm (Cho et al., 2014).

In this form, the model cell consists of three compartments: extracellular medium, cytoplasm and nucleus. We assume that the ratio of the cytoplasmatic volume to external volume (V _cyt /V _out ) is one, and that the radio of the cytoplasmatic volume to the nuclear volume (V _cyt /V _nuc ) is 2. In each compartment the concentration of the molecules are measured in pg mL⁻¹.

In the circuit of Figure 2, the rate of production of free cytoplasmic NF κB depends on the rate of its release from the IκB-NF κB complex, which is proportional to the product of the concentrations of IκB and IKK (Figure 3A), its own rate of production due to nuclear NF κB transcriptional activity inhibited by HDAC2 (Figures 1, 3A), the rate at which is recaptured by IκB (Figure 3A), and the rate of transport into the nucleus: d ( n f ) d t = V n f max p n f o n ( t ) + k a i k k ⋅ i κ b n f − k 4 n f ⋅ i κ b − D ⋅ n f (1) where nf is the concentration of free cytoplasmic NF κB, V n f max is the maximum rate of translation of the NF κB gene, p n f o n ( t ) is the probability that the NF κB gene is in its active state (on), ikk is the concentration of IKK, iκbnf is the concentration of IκB-NF κB complex, k _a and k ₂ are rate constants and D is the transport constant (Table 1).

We assume that the concentration of IKK activated by IL-6 is proportional to the external concentration of the cytokine: i k k = γ 1 ⋅ i l 6 , where γ ₁ is a constant of proportionality. In a similar form, PGE2 increases the rate of cleavage of the complex IκB-NF κB by phosphorylation of the IKKα subunit by Akt (Bai et al., 2009) (Figure 3B). Thus, we assume that the concentration of Akt is proportional to the external concentration of PGE2: a k t = γ 2 ⋅ p g e 2 . Substituting both expressions in Eq. 1 we finally obtain: d ( n f ) d t = V max n f p n f o n ( t ) + k 1 i l 6 ⋅ i κ b n f + k 11 p g e 2 ⋅ i κ b n f − k 4 n f ⋅ i κ b − D ⋅ n f (2) where k ₁ and k ₁₁ are the new rate constants (Table 1).

The rate at which the complex IκB-NF κB is broken in presence of IKK depends on the rate at which IKK and Akt remove IκB from the complex under the action of IL-6 and PGE2, and the rate at which the complex is formed again: d ( i κ b n f ) d t = k 4 n f ⋅ i κ b − k 5 i l 6 ⋅ i κ b n f − k 11 p g e 2 ⋅ i κ b n f (3) where iκbnf is the concentration of the complex IκB-NF κB, and k ₄ and k ₅ are rate constants (Table 1).

The rate at which NF κB enter the nucleus depends on the rate of transport of the molecule into de nucleus, adjusted for the change in volume between the cytoplasmic and nuclear compartments, and on the rate of inhibition of nuclear NK κB by its inhibitors: d n f ∗ d t = D ( V c y t V n u c ) n f − k deg n f ∗ (4) where nf* is the nuclear concentration of NF κB and k _deg is a rate constant.

The probabilities of activation of genes IL6, Cox2, IκB and NF κB (Figures 3A,B) are given by: d p I L 6 o n ( t ) d t = k 2 ( n f ∗ h d ∗ + 1 ) ( 1 − p I L 6 o n ( t ) ) − k − 2 p I L 6 o n ( t ) (5) d p I κ B o n ( t ) d t = k 3 ( n f ∗ h d ∗ + 1 ) ( 1 − p I κ B o n ( t ) ) − k − 3 p I κ B o n ( t ) (6) d p C o x 2 o n ( t ) d t = k 7 ( n f ∗ h d ∗ + 1 ) ( 1 − p C o x 2 o n ( t ) ) − k − 7 p C o x 2 o n ( t ) (7) d p N F κ B o n ( t ) d t = k 15 ( n f ∗ h d ∗ + 1 ) ⋅ ( 1 − p N F κ B o n ( t ) ) − k − 15 p N F κ B o n ( t ) (8)

Supplementary Material S2 for the mathematical deduction of these equations. For Eqs (5)–(7), the variable p i o n ( t ) , with i ∈ { IL - 6 ,  I κ B , Cox 2 , NF  κ B } , is the probability that the gene i is expressed at time t, hd* is the amount of HDCA2 in nucleus, and k ₂ , k _-2 , k ₃ , k _-3 , k ₇ , k_-7, k ₁₅ , k _-15 are rate constants (Table 1).

We assume that free cytoplasmic HDAC2 is produced according to the rate equation: d ( h d ) d t = r η ( n s p 5 + 1 ) − D 2 ⋅ h d (9) where hd is the amount of free cytoplasmic HDAC2 that is produced at a constant rate r, D ₂ is the rate of transport of free HDAC2 into the nucleus. In this equation, nsp5 is the amount of Nsp5 in the infected cell that sequester HDAC2 in cytoplasm (El Baba & Herbein, 2020), and η is a constant (Table 1). Free cytoplasmic HDCA2 enters the nucleus according to the equation: d ( h d ∗ ) d t = D 2 ( V c y t V n ) h d − k deg 2 h d ∗ (10) where hd* is the amount of HDAC2 in the nucleus, and k _deg2 is a rate constant (Table 1).

Nsp5 is a protein with only one link in the virus interactome (Gordon et al., 2020; Díaz. 2020a; Hekman et al., 2020), which is an input to the circuit of Figure 2. However, the real kinetic mechanism of production of this viral protein is unknown, thus we assume a plausible simple kinetic mechanism that leads Nsp5 to have a maximum constant concentration of 20 nM. We set the parameters ϕ and k ₁₈ to the values shown in Table 1: d ( n s p 5 ) d t = ϕ − k 18 n s p 5 (11) where ϕ is the constant rate of production of Nsp5, and k ₁₈ is a rate constant (Table 1). The effect of the variation of parameter ϕ on the dynamics of the circuit of Figure 2 is analyzed in Supplementary Material S3.

During the early stage of infection by SARS-CoV-2, the rate of variation of the concentration of IL-6 outside the lung cell depends on the rate of production of IL-6 by monocytes, leucocytes and other cells under stimulation by Tall-like Receptors (TLRs) [denoted by α (t)] (Jafarzadeha et al., 2020; Magro, 2020), on the rate of production by lung cells due to nuclear NF κB transcriptional activity enhanced by Nsp5 (Figure 3A), and on the rate of degradation of IL-6 by different mechanisms. Finally, we assume that the concentration of IKK activated by IL-6 is proportional to the external concentration of the cytokine (Eq. (2)): d ( i l 6 ) d t = α ( t ) + V I L 6 max p I L 6 o n ( t ) − k 5 i l 6 ⋅ i κ b n f − k 6 i l 6 (12) where V I L 6 max is the maximum rate of translation of IL6, p I L 6 o n ( t ) is the probability that the IL6 gene is activated at time t, and k ₅ and k ₆ are rate constants (Table 1).

The rate of variation of the amount of free IκB in the system is the balance between the rate of translation of gene IκB and the rate of formation of new IκB-NF κB complexes: d ( i κ b ) d t = V I κ B max p I κ B o n ( t ) − k 4 i κ b ⋅ n f (13) where V I κ B max is the maximum rate of translation of IκB, and p I κ B o n ( t ) is the probability that the gene IκB is activated at time t (Table 1).

The rate of variation of the amount of the enzyme Cox2 in the system depends on the rate of translation of the Cox2 gene, and its rate of inhibition by different mechanisms: d c o x 2 d t = V C o x 2 max p C o x 2 o n ( t ) − k 10 c o x 2 (14) where cox2 is the amount of the enzyme in the system, V C o x 2 max is the maximum rate of variation, p C o x 2 o n ( t ) is the probability that the gene Cox2 is activated at time t, k ₈ , k ₉ and k ₁₀ are rate constants (Table 1).

The rate of variation of the amount of PGE2 in the system depends on the amount of Cox2, on the amount of IL-6, and on its rate of degradation: d p g e 2 d t = k 17 c o x 2 − k 11 p g e 2 ⋅ i κ b n f − k 12 p g e 2 (15) where k ₁₁ , k ₁₂ and k ₁₇ are rate constants (Table 1).

In this section we present in detail and focus on the results that are relevant to understand the role of Nsp5 in the process of acute lung inflammation. The rest of the analyses and results related to the mathematical aspects of the model and its stability analysis are reported in Supplementary Material S3.

We solved the model using the Euler predictor-corrector and Runge-Kutta 4,5 methods with time step of 0.05 s and 720,000 integration steps to simulate a period of 10 h post infection, which is the estimated time for the first bursting of new virions into the extracellular medium (Bar-On et al., 2020; Kumar et al., 2020). We made a Boolean simulation of the circuit to confirm the dynamical behavior predicted by the ODEs model. We used the parameter values shown in Table 1 to generate an arbitrary reference state (ARS) of the system for the initial conditions iκbnf = 5 pg ml⁻¹, and the rest of the variables equal to 0 at time t = 0. In this case, the input α (t) = 10 pg ml⁻¹s⁻¹ and nsp5 = 0 for t ∈ [ 0 , ∞ ) . In this ARS, the translation velocities of all genes are set to 2 pg mL⁻¹ s⁻¹ (Hausser et al., 2019). In the case of the ODEs model, both numerical methods produce identical outputs indicating no bias in the predicted dynamics of the circuit due to the integration method and that the transitory sharp initial increase observed in the concentration of the proteins of the circuit are to due to the transitory initial response of the feedback loops to the input. The solution of the Boolean model confirms that the trajectories of the circuit lead to a stable node or attractor where all the nodes are in on state, except HDAC2 that is in state off, when Nsp5 is in state on.

Figure 4A shows that the output of the model is the transitory activation of NF κB, Cox2 and PGE2, in contrast with a sustained high activation of IL-6. In order to clarify the source of the steady activation of IL-6, we solved the model with the same initial conditions and α (t) = 0. In this case, the system does not turn on during all the time of simulation (data not shown), indicating that the signal from monocytes and other cells is a necessary condition to start the IL-6 sustained production in lung cells (Figure 4A). Figure 4B shows that when nsp5 ≠ 0, α(t) = 0 and equal gene translation rates, the response of the system is a high sustained level of activation of NF κB, Cox2, and PGE2 with respect to the ARS. In this case, the concentration of IL-6 in the cell was only slightly affected by the virus. Figures from 4A to 4D show the first 5 minutes of simulation during which the circuit of Figure 2 reaches its steady state. This result indicates that the presence of the virus induce a sustained and persistent production of NF κB, Cox2, and PGE2 in the circuit of Figure 2.

Translation rates of the set of genes of the model rarely have the same value, even when their activation is promoted by the same transcription factor (NF κB in this case). As shown in Table 1, the results of Figures 4A,B are obtained when the rates of translation are the same for all genes. However, in real cells this is not the case, and Cox2, Il-6, NF κB and IκB are translated at different rates. In Homo sapiens the estimated average rate of translation ( V max ) of a gene is about 10⁴ proteins per hour (Hausser et al., 2019). However, this average value can be increased according to the state of activity of the cell. In the present work, we arbitrary chose the value of V max = 2  pg mL − 1 s − 1 for all genes in the ARF (Table 1). However, we made parameter variation to know how different values of this set of parameters affect the qualitative behavior of the model (Supplementary Material S1). Figure 4C shows that when in the system V n f max = 10  pg mL − 1 s − 1 , V I κ B max = 15  pg mL − 1 s − 1 , V C o x 2 max = 2  pg mL − 1 s − 1 , and V I L 6 max = 2  pg mL − 1 s − 1 , α(t) = 10 pg ml⁻¹s⁻¹ and nsp5 = 0, the circuit tends to a steady state in which NF κB, Cox2, and PGE2 exhibit a constant low concentration during the time of simulation, while IL-6 steady concentration is unaffected with respect to the value shown in Figure 4A. When nsp5 ≠ 0, the steady concentration of NF κB, and Cox2, are increased ∼10 fold with respect to their maximum concentration in the ARS, while PGE2 concentration is increased ∼4 fold; Nsp5 produces a ∼4 fold decrease in IL-6 steady concentration with respect to the reference concentration shown in Figure 4A.

This result indicates that a the difference in translation rates between NF κB and IκB (the intensity of the inhibition of NF κB by IκB, Figure 3A) in absence of the virus with respect to their values in ARS induces a low sustained production of NF κB, Cox2, and PGE2 without affecting IL-6 production. However the presence of Nsp5 induces a high production of NF κB and Cox2, and a moderate increase in PGE2 concentration that decreases IL-6 concentration, indicating that a higher rate of IκB translation with respect to NF κB translation can decrease the concentration of the cytokine IL-6 in presence of the virus probably due to a lower concentration of PGE2in the circuit (Figures 2, 3A).

Figure 4D shows that when α (t) = 10 pg ml⁻¹ s⁻¹, nsp5 = 0, V n f max = 10  pg mL − 1 s − 1 , V I κ B max = 15  pg mL − 1 s − 1 , V C o x 2 max = 2  pg mL − 1 s − 1 , and V I L 6 max = 20  pg mL − 1 s − 1 the circuit tends to a steady state in which NF κB, Cox2 and PGE2 exhibit a constant low concentration during the time of simulation, while IL-6 is unaffected with respect to its steady concentration in Figure 4A. However, when nsp5 ≠ 0 the dynamical behavior of the circuit changes because NF κB and Cox2 show the same steady concentration that in Figure 4C, but PGE2 has a higher steady concentration. In this case, IL-6 shows an increase of ∼2 fold with respect to the level shown in Figure 4C. Thus, a high rate of IL-6 translation can overcome the effect of an increased rate of IκB translation in presence of the virus, allowing an overproduction of the cytokine.

All these results suggest an important dynamical property of the circuit shown in Figure 2: a high value of the parameter V I L 6 max produces an over production of NF κB, Cox2, PGE2, and IL-6 in the circuit in presence of Nsp5 with respect to the ARS (Figures 4A,D). We name over stimulated immune state (OSIS) to the state of the system shown in Figure 4D, for which: V n f max = 10 pg ml⁻¹ s⁻¹, V I κ B max = 15 pg ml⁻¹ s⁻¹, V C o x 2 max = 2 pg ml⁻¹ s⁻¹, V I L 6 max = 20 pg ml⁻¹ s⁻¹, nsp5 ≠ 0 and α (t) = 10 pg ml⁻¹ s⁻¹.

Figure 5 shows the response of the circuit to a step function defined as: α ( t ) = { 10 0 ≤ t < 10 0 t > 10 (16)

The latter simulates an external signal of IL-6 initially generated at a rate of 10 pg ml⁻¹ s⁻¹ during 10 min, and switched off after the 10 min. Figure 5A is the response of the system in ARS to the step signal, in absence and presence of Nsp5. Interestingly, in the absence of Nsp5, the circuit exhibits a low concentration of NF κB, Cox2 and PGE2, and a high concentration of IL-6, all which are zero when the external signal is switched off. In presence of Nsp5 the ARS is perturbed by α (t), and reaches a new steady state with high concentrations of NF κB, Cox2, and PGE2 that remain even when the external signal is switched off. Furthermore, the high IL-6 concentration reached when the circuit is turned on falls to a lower concentration different from zero when α (t) = 0.

In order to reduce the effect of Nsp5 on the circuit of Figure 2, we tested the effect of Nimesulide on NF κB production. Nimesulide, a Cox2 inhibitor, acts inhibiting the synthesis of Cox2 and lipooxigenase enzyme and their products (Suleyman et al., 2008). We modified Eq. 14 to introduce an inhibition term that decreases the rate of production of Cox2 in presence of Nimesulide: d c o x 2 d t = V C o x 2 max p C o x 2 o n ( t ) n i m + 1 − k 10 c o x 2 (17) where nim is the amount of Nimesulide at time t, whose rate of variation of its concentration in the cell is: d ( n i m ) d t = r 2 − k 19 n i m (18) where r ₂ = 40 pg ml⁻¹ s⁻¹ and k ₁₉ = 3 s⁻¹. Solution of the modified model (Eq. (17)) is presented in Figure 6. In Figure 6A, Nimesulide is applied to the circuit in presence of Nsp5, and with all gene translation rates set to 2 pg ml⁻¹ s⁻¹. In this case, Nimesulide has no effect on NF κB, a slightly effect on IL-6, and a significant decrease in Cox2 and PGE2 concentrations to values near to zero in presence of 14 pg ml⁻¹ of the drug. The same effect is observed when the circuit is in OSIS (Figure 6B). Thus, Nimesulide is not an effective anti-inflammatory agent in the context of the dynamic circuit modeled here, and does not impact the concentration of IL-6 in the presence of Nsp5, even at high concentrations.

Dexamethasone (DX), is a potent synthetic glucocorticoid that inhibits NF κB transcription, and is another anti-inflammatory drug that is being used to ameliorate acute inflammatory states due to SARS-CoV-2 infection. In order to test the effect of DX on the inflammatory process in the circuit of Figure 2, we modified Eq. 8 to introduce an inhibition term that decreases the probability of activation of NF κB: d p N F κ B o n ( t ) d t = k 15 ( n f ∗ ( h d ∗ + 1 ) ( D X + 1 ) ) ( 1 − p N F κ B o n ( t ) ) − k − 15 p N F κ B o n ( t ) (19) where DX is the amount of DX at time t, whose rate of variation of its concentration in the cell is: d ( D X ) d t = r 3 − k 20 D X (20) where r ₃ = 30 pg ml⁻¹ s⁻¹ and k ₂₀ = 3 s⁻¹.

Figure 7 shows the effect of DX on the circuit. Figure 7A shows that when DX is applied to the circuit in presence of Nsp5, and when all genes have the same translation rate of 2 pg ml⁻¹ s⁻¹, the drug eliminates NF κB, Cox2 and PGE2 from the circuit, but has a null effect on the concentration of IL-6. However, when the circuit is in OSIS (Figure 7B), DX has little effect on the immune response in presence of Nsp5. Figure 7C shows that when all translation rates are the same (2 pg ml⁻¹ s⁻¹), DX can shut down the immune response in presence of the virus, when the circuit input α (t) is the function of Eq. 16. However, DX cannot shut down the effects of Nsp5 when the circuit is in OSIS, having a little anti-inflammatory effect (Figure 7D). Hence, our results suggest that DX is an effective anti-inflammatory drug during mild inflammatory states, but is less effective during critical acute lung inflammation due to SARS-CoV-2 (Horby et al., 2021). This dynamical behavior of the circuit of Figure 2 is unaffected by higher doses of DX. When DX and Nimesulide are applied together (with 35 pg ml⁻¹ of DX and 40 pg ml⁻¹ of Nimesulide) to the circuit in OSIS, in which α (t) is the function of Eq. 16, there is also a limited anti-inflammatory effect (Figure 8A).

In contrast, the predicted limited effect of Nimesulide and DX on reverting the OSIS, the molecular docking of Nsp5 with a hypothetical drug like Saquinavir (Xu et al., 2020), is able to completely eliminate inflammation. We tested this hypothesis modifying Nsp5 production in Eq. (11) with an inhibitory term: d ( n s p 5 ) d t = ϕ d r u g + 1 − k 18 n s p 5 (21) where drug is the amount of the hypothetical drug that can remove Nsp5 from the circuit, whose concentration at time t is given by: d ( d r u g ) d t = r 4 − k 21 d r u g (22) with r ₄ = 30 pg ml⁻¹ s⁻¹ and k ₂₁ = 3 s⁻¹.

Figure 8B shows that when the hypothetical drug is applied to the circuit in OSIS, with α (t) given by Eq. 16, the inflammatory process is completely shut down, suggesting that Nsp5 is effectively the key cause of the acute lung inflammation during SARS-CoV-2 infection, and a fundamental potential target for its treatment.

The complexity of the immune response to infections implies a complex set of different kinds of cells that involve complex regulatory networks and molecules. Together, they underlie specific responses that may lead to resolution or to inflammatory responses whose dynamical properties depend on the networks structure (i.e., the number of nodes, links, and the distribution of links between the nodes) occurring at different scales and levels (Lipniacki et al., 2004; Morel et al., 2006; Eftimie et al., 2016). The structure of such complex and dynamical network is not random due to the specific nature of its nodes and links, so that if the network has a hierarchical structure its hubs determine the overall network dynamics and the structure of its phase space (Cessac, 2010; Eftimie et al., 2016; Díaz. 2020b). In high dimensional phase spaces, steady points of diverse kinds can coexist with isolated closed curves (limit cycles), strange attractors of fractal structure, and other limit sets giving rise to a variety of biological dynamical behaviors like molecular switches, periodic and quasi-periodic oscillations, and bursting (among others). Two central problems in the analysis of this kind of dynamical systems are: 1) the structure of the phase space is strongly dependent of the value of the set of parameters of the dynamical system, and 2) the nonlinearity characteristic of biological systems such as the one treated here (Supplementary Material S1) (Nayak et al., 2018).

In most of the ODEs-based biological models so far studied, there is limited knowledge of the real values of the parameters, and this leads to uncertainty concerning the extent to which the model explains the real process under study. In this case, the validation of the model relied on experimental data. Nonetheless, with the exception of some particular cases (Hodgkin and Huxley, 1952), qualitative models have to be proposed and quantitative validations remain ahead (Eftimie et al., 2016). This challenge increases with the dimensionality of the model. Thus, is necessary to choose some particular properties of the network under study (inflammation network in this case) that incorporate key dynamical features of the whole system under study, and that can also aid to reduce the number of parameters of the model to the minimum possible.

In this sense, we limited the scope of this work to the study of the inflammation process of epithelial lung cells during SARS-CoV-2 infection. The epithelial lung cells form a sub-network that have two highly connected nodes: the internal signal integrator NF κB, and the external signal integrator IL-6 Receptor (IL6R) (Magro, 2020), which is activated by external IL-6 (Figures 2, 3). Both integrators trigger the activation of two important intermediate processes: Cox2 and PGE2 synthesis. Thus, we propose, as a first approximation, that the dynamical behavior of this these four nodes of the circuit of Figure 2 reflect with certain accuracy the dynamical properties of the complex inflammation sub-network of lung cells. Thus, the dynamics of IL-6, NF κB, Cox2 and PGE2 nodes is described by a model that consist of 12 nonlinear coupled ordinary differential equations (Table 1), that we assume represents the central dynamical core involved in the lung cells inflammation process in response to SARS-CoV-2 infection. Such core includes three coupled feedback loops shown in Figure 3. We used this low-dimensionality model to explore the effect of Nsp5 in the qualitative dynamics of the inflammation circuit of Figure 2 (Díaz, 2020b; Kumar et al., 2020).

SARS-CoV-2 main protease Nsp5 has a role as an epigenetic regulator of host DNA expression in lung cells; and thus its physical association with HADC2 (Gordon et al., 2020; Kumar et al., 2020) avoids the movement of the deacetylase from the cytoplasm to the nucleus, increasing the probability of the binding of NF κB to the promoter site of its target genes (Eqs. 5–8) (El Baba and Herbein, 2020). NF κB is a central integrator of the signals that initiate the inflammation process (Figure 1), which include IL-6 secreted by monocytes, leucocytes, macrophages, among other cells, together with signals from TLR4 and TLR7 that suppress the IκB-NF κB complex in a MyD88 independent manner (Vidya et al., 2017). NF κB enhanced transcriptional activity includes the transcription of IL-6, Cox2, IκB and NF κB itself that are key genes in the process of inflammation (Figure 3).

Results obtained from our model suggest that Nsp5 can effectively transform a weak inflammation response (Long et al., 2020) into a persistent sustained one that has a relative high concentration of the cytokine IL-6 (Figure 4B). In the example of Figure 4B, all genes have the same translation rate, a condition that have low probably of occurrence in the real human immune system in which every gene is subject to different processes of post transcriptional regulation (Hausser et al., 2019). When this restriction is removed, and IL-6 has a high translation rate, Nsp5 boosts the concentration of the four main proteins of the circuit with respect to their respective values in ARS, and the inflammation response becomes strongly persistent with a high IL-6 concentration. Furthermore, in this state that we named OSIS, the dynamics of the circuit of Figure 2 becomes independent of any kind of external signal and never turns off, leading to the deregulation of the production of the cytokine IL-6 (Figure 5B). In this form, Nsp5 changes the qualitative dynamical behavior of the system, in which the trajectories that initially span around the ARS stable node now span around the OSIS stable node (Supplementary Material S3) (Díaz, 2020b).

An important biological consequence of this change in the qualitative dynamics of the circuit is that ACE2 receptor synthesis also becomes persistent and autonomous of any external signal to the lung cells (Figure 2). Furthermore, neither Nimesulide nor Dexamethasone or both can completely eliminate the OSIS suggesting that this dynamical behavior could be the form in which the virus assures the persistence of its reproductive cycle without any perturbation from the natural defenses of the body. This could be also a possible clue about the possible adaptive substitution of N protein used in SARS-CoV for the use of the main protease Nsp5 in SARS-CoV-2 as the switch for the inflammation process, which is necessary for ACE2 sustained production, taking into consideration that SARS-CoV N protein targets Cox2 while SARS-CoV-2 Nsp5 protein targets the master integrator NF κB. If this hypothesis is true, this is the probable cause of the increased pathogenicity of SARS-CoV-2 with respect to SARS-CoV.

It is of interest that the persistence of the OSIS is supported by a high translational rate of IL-6, which in the model must have a value of V I L 6 max ≥ 7  pg mL⁻¹s⁻¹. This result pictures a probable hypothetical scenario where the external signal α (t) triggers the initial inflammation response with a low translation rate of IL-6 that can be shut down with DX. As viral infection continues, IL-6, NF κB, and IκB increase their rate of translation ( V C o x 2 max value has a little weight in this process) and the dynamical behavior of the circuit of Figure 2 becomes independent of the value of α (t), and self-sustained, i.e., the control of the production of the cytokine IL-6 is translated from the external immune cells to the lung cells, and cannot be completely stopped even with high doses of DX. In this form, in the model, V I L 6 max is a bifurcation parameter that changes the qualitative dynamical behavior of the circuit of Figure 2 when V n f max = 10  pg mL⁻¹s⁻¹ and V I κ B max = 15  pg mL⁻¹s⁻¹ are constant values (Supplementary Material S3).

In this extreme situation (OSIS), the unique treatment possible suggested by our model is an inhibitor of Nsp5 like Saquinavir (Xun et al., 2020). In this case, the OSIS is completely shut down, which indirectly implies the down production of ACE2. Saquinavir is an anti-retroviral drug used against the main protease of HIV-1 (Bensussen et al., 2018; Xun et al., 2020), with some undesirable effects like diarrhea, abdominal pain, and nausea. We suggest that this drug could be a plausible treatment against the effects of SARS-CoV-2 in acute lung inflammation, and as a complement to the new vaccine against this coronavirus.

Interferon (IFN) I and III are key secretory factors in the fight against viral infections in lung cells. During the early response to SARS-CoV-2 infection, interferon pathway is blocked by viral proteins like Nsp14, Nsp16, Nsp15, Orf9b, and Nsp3, among others, that interfere with MAVS and TLRs sensing, STATS 1 and 2 activation and IFN synthesis. Experimental results show that this downregulation of IFN activity correlates with an increased severity of COVID19 symptoms (Kim and Shin, 2021). Therefore, we did not include IFN pathway in our model as a primary cause of lung inflammation. However, our model suggest that the initial downregulation of IFN pathway could be partially counteract by the over activation of NF κB caused by Nsp5. The master regulator NF κB targets IFNs genes (Pfeffer, 2011; Kim and Shin, 2021) probably producing a late sustained production of IFN after the initial viral infection adding a component to the cytokine storm.

In a previous work (Díaz, 2020a) the role of the viral protein Orf8 as a hub of the SARS-CoV-2 infection was suggested. This highly connected protein has as main targets those processes related to the vesicular trafficking required for the ensemble of new virions (Figure 1) (Gordon et al., 2020). Thus, the inhibition of the effects of Orf8 in the host cells could block the reproduction cycle of the virus but not the inflammation process because Nsp5 is not directly linked to Orf8 (Gordon et al., 2020). Rapamycin has been suggested as a drug against the effects of Orf8 but produces severe immunosuppressant effects, which, according to the results of our model, will not be of care in cases of severe lung inflammation because Nsp5 uncouple the circuit of Figure 2 from the process of viral replication (Figure 1 and Figure 5B). In this case, a treatment with Rapamycin and Saquinavir can be an alternative for patients with severe lung inflammation. This suggestion needs experimental verification and clinical evaluation.

In this work we propose a dynamical model that puts forward a novel hypothetical mechanisms underlying lung cells inflammation process in response to SARS-CoV-2. In this scenario the main protease Nsp5 effectively enhances the inflammatory process, increasing the levels of NF κB, IL-6, Cox2, and PGE2 with respect to the reference state. When the translation rates of NF κB and IκB are increased to a high constant value, and the translation rate of IL-6 is increased above the threshold value of 7 pg ml⁻¹s⁻¹ the circuit enters in a persistent over stimulated immune state (OSIS) with high levels of the cytokine IL-6. The OSIS never shuts down by itself, and becomes autonomous of the signals from other immune cells like macrophages and lymphocytes. Our model explains why DX or Nimesulide have little effect on the OSIS, and the only means to suppress such acute and sustained inflammation state is by inhibiting the protein Nsp5 with drugs such as Saquinavir.

In this form, the model suggests, in accordance to our hypothesis, that Nsp5 is effectively a key node underlying severe acute lung inflammation during SARS-CoV-2 infection. The persistent production of IL-6 by lung cells during infection can be one of the causes of the cytokine storm observed in critical patients with COVID19. From an evolutionary point of view, the use of Nsp5 as the switch to start inflammation, and the consequent overproduction of the ACE2 receptor, could have driven the increased pathogenicity of SARS-CoV-2 with respect to SARS-CoV.

In the absence of quantitative data on the kinetics of Nsp5 in the host cell that could give us a better overview of the dynamical behaviour of the protein in the model proposed in this work, we highlight the importance of future in vitro and in vivo experimental research. In order to experimentally validate the predictions of our model is necessary to measure the rate of production of Nsp5 during the infection of the host lung cells, and compare it with the rate of production of other viral proteins with immunomodulatory effects as Orf 7b and Orf 8, to know the real magnitude of its contribution to the inflammatory process. Furthermore, the experimental inhibition of Nsp5 could be of aid to eliminate the inhibition of HDAC2, either catalytic or mimetic, in infected cells to confirm experimentally if the predictions of the model are in concordance to the real mechanisms that contribute to the cytokines storm. The implementation of techniques as plasmid induction or transfection for the unique expression of Nsp5 and its effectors in lung epithelial tissue also promises to be a simple option to validate the proposed mathematical model, and in turn it would be an ideal platform to measure the efficacy of Nsp5 inhibitors proposed by bioinformatics and clinics.

The present theoretical work is part of the project “Dynamics of the SARS-CoV-2 network” and is based on the experimental data reported in the literature at the moment.