key: cord-0845229-fd0zcfg5
authors: Romano, Alessandra; Casazza, Marco; Gonella, Francesco
title: Addressing Non-linear System Dynamics of Single-Strand RNA Virus–Host Interaction
date: 2021-01-15
journal: Front Microbiol
DOI: 10.3389/fmicb.2020.600254
sha: 74bc3733d6248a3d45ce93394413f48677b89713
doc_id: 845229
cord_uid: fd0zcfg5

Positive single-strand ribonucleic acid [(+)ssRNA] viruses can cause multiple outbreaks, for which comprehensive tailored therapeutic strategies are still missing. Virus and host cell dynamics are tightly connected, generating a complex dynamics that conveys in virion assembly to ensure virus spread in the body. Starting from the knowledge of relevant processes in (+ss)RNA virus replication, transcription, translation, virions budding and shedding, and their respective energy costs, we built up a systems thinking (ST)–based diagram of the virus–host interaction, comprehensive of stocks, flows, and processes as well-described in literature. In ST approach, stocks and flows are expressed by a proxy of the energy embedded and transmitted, respectively, whereas processes are referred to the energy required for the system functioning. In this perspective, healthiness is just a particular configuration, in which stocks relevant for the system (equivalent but not limited to proteins, RNA, DNA, and all metabolites required for the survival) are constant, and the system behavior is stationary. At time of infection, the presence of additional stocks (e.g., viral protein and RNA and all metabolites required for virion assembly and spread) confers a complex network of feedbacks leading to new configurations, which can evolve to maximize the virions stock, thus changing the system structure, output, and purpose. The dynamic trajectories will evolve to achieve a new stationary status, a phenomenon described in microbiology as integration and symbiosis when the system is resilient enough to the changes, or the system may stop functioning and die. Application of external driving forces, acting on processes, can affect the dynamic trajectories adding a further degree of complexity, which can be captured by ST approach, used to address these new configurations. Investigation of system configurations in response to external driving forces acting is developed by computational analysis based on ST diagrams, with the aim at designing novel therapeutic approaches.

Positive single-stranded ribonucleic acid [(+)ssRNA] viruses, including picornaviruses, flaviviruses, Togaviridae, and human coronaviruses (CoVs) (Ahlquist et al., 2003; Lam et al., 2016; Scutigliani and Kikkert, 2017; Primadharsini et al., 2019) , cause multiple outbreaks, for which tailored antiviral strategies are still missing (Zumla et al., 2016; Dinesh et al., 2020; Gordon et al., 2020b) . (+)ssRNA viruses package their genomes as messenger sense, single-stranded RNA and replicate those genomes solely through RNA intermediates in the cytosol of the host cells (Den Boon et al., 2010) . RNA-dependent RNA polymerases lack coreplicative and postreplicative fidelity-enhancing pathways; this final RNA genome copies incorporate mutations at a high rate (Lauring et al., 2013; Acevedo et al., 2014) , providing the viral quasi-species with a higher probability to evolve and adapt to new environments and challenges during infection (Burch and Chao, 2000; Vignuzzi et al., 2006) . The diversity is essential for both viral fitness (Wargo and Kurath, 2012) and pathogenesis because of the complex relationships among virus replication (VR), host cells, and immune system, as almost all (+)ssRNA viruses can delay antiviral innate immune response (Kobasa et al., 2007) in multiple ways (Diebold et al., 2003; Hogan et al., 2004; Meylan et al., 2005; Saito et al., 2008; Yang et al., 2015; Beachboard and Horner, 2016; Nelemans and Kikkert, 2019) . Host immunogenetic factors can be sensitive to a variation in the viral load, leading to a defective response of the innate immunity, that could explain the variable clinical course of infection (Fanning et al., 2001; Nelemans and Kikkert, 2019) .

Recent studies confirmed the complexity of viral dynamics, whose fitness is improved by the complex interactions with the host proteins (Bosl et al., 2019; Sruthi and Prakash, 2019) , as previously described in modeling the virus-host interactions at subcell and cell levels (Dapat and Oshitani, 2016; Jonsdottir and Dijkman, 2016; Gao et al., 2017; Patzina et al., 2017; Gordon et al., 2020b) . However, models that address a specific aspect of the virus-host interaction do not capture the wide range of intertwined spatial and temporal (hours to days) dynamic scales (Apweiler et al., 2018) , which are related to the interactions of different concurrent hierarchical levels. For this reason, we aimed at describing the virus-host cell interaction as a dynamic system by a systems thinking (ST) approach (Northridge and Metcalf, 2016) .

In ST, the behavior of a dynamic system can be described and predicted by the temporal evolution of its configurations, given by hierarchical feedback loops and self-organization. The configuration evolution then can be analytically computed by proper simulators (Odum and Odum, 2000; Tegner et al., 2009; Hassmiller Lich et al., 2017; Spill et al., 2018) to further address suitable leverage points for intervening (Meadows and Wright, 2008) . In particular, dynamic models, where the temporal evolution of extensive variables (stocks) is simulated in the form of trajectories, derive their initial conditions from available and assessed evidence. Varying the key system parameters, trajectories represent the possible evolutive (structural) patterns of the system at issue, becoming abstracted with respect to local specific attributes related to single case studies. However, being stocks associated to system observables, the relation with those attributes is maintained, making it possible to compare predicted trajectories with observed data and constituting suitable counterfactuals with respect to the laboratory measurements results.

System dynamics (SD) approach is mostly used for strategic modeling, typically for ecological and socioeconomic systems, to understand the supply chain performance. In particular, results of SD modeling provide a set of alternative evolutive patterns in form of graphs, capturing the internal dynamics of a system even in lack of some experimental data to fit. These can provide alternative scenarios that, when fitting experimental evidences, indicate the most effective leverage points to control the system evolution.

In this article, we show that, by approaching the host-virus interaction as a dynamic systemic problem (Sterman, 2002; Meadows and Wright, 2008) , it is possible to identify potential systemic leverage points to minimize the release of virions, so addressing effective systemic intervention strategies.

The basic SD element is the stock and flow diagram. Stocks are countable extensive variables Q i , i = 1, 2,. . ., n, relevant to the study at issue, that constitute an n-ple of numbers (possibly derived from experimental measurements), which at any time represents the state of the system. A stock may change its value only upon its inflows and/or its outflows, represented by arrows entering or exiting the stock. Processes are any occurrence capable to alter-either quantitatively or qualitatively-a flow, by the action of one or more of the system elements. In a stationary state of the system, stock values are either constant or regularly oscillating. Processes, which cause the stationarity or perturbation of a system, must be activated by a driver, acting on the flows where the process is located. The pattern of the feedbacks acting in the system configurations is the feature that ultimately defines the systems dynamics. Each flow depends on the state variables Q i by relationships of the kind dQ i /dt = kf (Q j ), i, j = 1,. . ., n, where n is the number of stocks in the system.

The stocks and flows inventory reported in Table 1 was based on information from existing knowledge on the biological mechanisms at issue, listing the variables and parameters necessary to set up the equations describing the system dynamics. Turnover times of stocks included in the RNA-virus-host interaction ST diagram have been reported in Table 2 , derived from the available literature (Konig et al., 2008; Friedel and Haas, 2011; Pfefferle et al., 2011; Jourdan et al., 2012; Munday et al., 2012; Naji et al., 2012; Wu et al., 2012; Emmott et al., 2013; Garcia-Dorival et al., 2014; Verchot, 2014; Watanabe et al., 2014; York et al., 2014; Zheng et al., 2014; Dong et al., 2016; Kuo et al., 2016 Kuo et al., , 2018 Wang et al., 2016 Wang et al., , 2017a Gao et al., 2017; Hafirassou et al., 2017; Khamina et al., 2017; King et al., 2017; Martinez-Gil et al., 2017; Patzina et al., 2017; Coyaud et al., 2018; Iwasaki et al., 2018; Lescar et al., 2018; Ziegler et al., 2018; Bosl et al., 2019; Chakravorty et al., 2019; Garcia-Moreno et al., 2019;  Odum, 1996 Odum, , 2002 Odum and Odum, 2000 Wheatley et al., 1980; Eden et al., 2011; Siwiak and Zielenkiewicz, 2013 Q 2B Long-half-life proteins Wheatley et al., 1980; Eden et al., 2011; Siwiak and Zielenkiewicz, 2013 J 3 + J 4 + J 5 + J 6 + J 7 + J 50 vRNA replication and virion shedding, virus production after infection 6 h Sedmak and Grossberg, 1973; Baccam et al., 2006 Rothenburg and Brennan, 2020). All stocks, flows, and processes were expressed using a common proxy unit, representing the energy embedded, transmitted, and used, respectively, during the system operation. The proxy unit was expressed as the number of ATP (and ATP-equivalent) hydrolysis events (Mahmoudabadi et al., 2017) . This choice allowed calculating each parameter of the system on the basis of stocks and characteristic times of the flows derived from the literature, without further need for experimental data.

After setting the initial conditions at time 0 for the stocks, system solutions were obtained using recursive computation for a relative short period of time (identified with the median life of an epithelial cell, 7 days), in order to appreciate the model dynamic behavior. The computational model based on a set of differential equations that describe the rates of change of all stocks in the ST diagram (Odum and Odum, 2000; Bossel, 2007) was developed using the open-source software package SCILAB 1 , which uses approximation techniques to evaluate stocks. Given a set of initial conditions for the stocks (i.e., the initial state of the system) and a set of phenomenological coefficients k associated to flows, the set of interconnected equations was treated by a standard finite-different method, taking care of choosing a time step short enough to evidence the dynamics of any of the studied processes. The coefficients k i were calculated, on the basis of literature data, considering the dynamics of each single stock, by quantifying flows and stocks during the time interval set as simulation step, as shown in Table 1 . When different flows coparticipate in a process, each coefficient gathers all the actions that concur to the intensity of the outcoming flow(s). In detail, the parameters used to run the model (i.e., the set of values for the k i coefficients), describing the reaction of each system component to a change in any other one, were derived from the stocks turnover times (Odum and Odum, 2000; Bossel, 2007) . Therefore, the host-virus interaction computational model, built on experimental evidences as listed in Tables 1, 2, is not specific for a unique virus, but may represent the patterns of any virus-host interaction, in which stocks, flows, and processes are those relevant for the operation of the system at issue.

The reliability of both available data and modeling was tested by evaluating the effect of the variation of each of the most relevant input data (stocks and processes) on the system trajectories. Unfortunately, here is not a single comprehensive procedure suitable for the validation of all dynamic models, being dependent on their usefulness, in turn referred to the very purpose of the model itself (Grüne-Yanoff and Weirich, 2010). We chose the sensitivity analysis approach (Qudrat-Ullah, 2012; Hekimoglu and Barlas, 2016) , which allows to see to what extent a variation on these values can lead to alternative evaluations of the system dynamics. In particular, we applied a 50% variation (either positive or negative) to those parameters that the results were more sensitive to. As expected, while the corresponding simulations varied as well, the general patterns presented in the following remained the same, especially concerning the overall trends shown by comparing the groups of simulations, providing a model validation.

First, we identified the important structures in the system and then used to build up the stock-flow diagram of the virus-host interaction system. In Figure 1 , symbols were borrowed from the energy language (Odum and Odum, 2000; Brown, 2004) : shields indicate the stocks; big solid arrows indicate the processes; line arrows indicate the flows; dashed lines show the controls exerted by the stocks on the processes.

The dynamics of energy allocation for protein synthesis contained in the stock Q 1 depended on the cell bioenergetics, e.g., the number of mitochondria, OX-PHOS activity levels, and cell cycle phase (Murayama et al., 2008; Canto et al., 2009; Li et al., 2020) . In the absence of virus, the stationary configuration was given by energy required to flow from stock Q 1 (via J 1 , J 2A , and J 2B flows) to stocks Q 2A and Q 2B to produce, respectively, short-and long-half-life proteins, which could, in turn, be recruited by VR machinery. J 20A and J 20B , grouped into the flow J 20 , represented the outflow of folded, fully functional proteins addressed to secretion or surface exposure. Based on basal proteostasis of host cell, recovery of energetic sources from proteins not addressed to leave the system could be possible via several complex processes (e.g., proteasomal degradation, and autophagy), identified by flows of materials J 21A and J 21B , respectively, from Q 2A and Q 2B back to Q 1 .

The viral load in the system, expressed by the stocks Q 3 (identified as viral RNA content to be used for viral transcription and translation), Q 4 (translated viral proteins content), and Q 5 (full assembled virions to shed virus outside), diverted, at the time of infection, resources directly from Q 1 (through flows J 13 , J 15 , and J 17 ) and Q 2B (through flows J 23 , J 25 , and J 27 ). Virions shedding was represented by the flows J 7 and J 50 through the contribution of the host flows J 17 and J 27 . The output flows J 4 and J 5 were set to be effective only if the value of the respective stock Q 3 and Q 5 was higher than a threshold, as represented by the two switch symbols in the diagram.

We identified four feedback loops (represented by dot lines in Figure 1 ): (i) the positive control of Q 2A stock on the energy supply process (occurring when more structural host proteins operate to maintain the energetics homeostasis of the host cell); (ii) the positive control of Q 3 stock on the VR process (highlighting that the more viral RNA is in the system, the more intensive replication can occur if host sources are available); (iii) the positive control of Q 4 stock on the processes of synthesis and maturation of host proteins (highlighting that the more viral proteins are made, the more host proteins are synthesized to be recruited in the virion assembly machinery, increasing J 2B ); (iv) the positive control of Q 5 stock on the VR process (highlighting that the more virions are produced, the more resources are diverted from the host cell to viral replication).

First, a computational model was derived from the stock-flow diagram shown in Figure 1 using the standardized workflow of systemic modeling (Odum and Odum, 2000; Xue et al., 2018) . Figure 2 shows two different system self-organized patterns (configurations) to guide reader in the overall comprehension of the proposed approach. The virus-host interaction was represented as an evolving set of simulated trajectories, to which the positive value of Q 3 stock had given access, using preexisting stocks, processes, and flows of the host cell, followed overtime by progressive filling of Q 4 and Q 5 stocks. In (Figure 2A ) configuration, the viral load is null (the stocks Q 3 , Q 4 , and Q 5 are empty), and the values of stocks Q 1 , Q 2A , and Q 2B are constant; thus, the system behavior is stationary (Figure 2B) . At time of infection, the Q 3 stock was fed, and its proteins could interact with the host proteome to sustain RNA replication. Based on previous works in the field (Wei et al., 1995; Adelman et al., 2002; Mohler et al., 2005; Regoes et al., 2005; De Boer et al., 2010) , we identified a time delay of 2-6 h required to record changes in the Q 5 stock.

Moreover, the value of Q 5 varied over time due to changes that occurred at different timepoints in the stocks Q 2B , Q 3 , and Q 4 . Thus, the network of flows and feedbacks could identify a new configuration (Figure 2C ), to generate a non-stationary FIGURE 1 | The energy systemic diagram of a cell infected by ss+ RNA virus. Stock-flow diagram of the virus-host interaction system. In the upper right box are the meaning of symbols. The color code is as follows: blue for host cell energy stocks and relative inflows and outflows; red for virus energy stocks and relative inflows and outflows; green for external energy inputs and external driving forces F corresponding to different therapeutic strategies. The lower box lists the biological contents of the stocks, all expressed in terms of energy (ATP-equivalent units).

behavior ( Figure 2D) , where the values of Q 3 , Q 4 , and Q 5 stocks evolved in a non-linear way (Supplementary Figure 1) , to maximize the value of virions stock in the configuration ( Figure 2E) . We define the (Figure 2A ) configuration as healthy, the ( Figure 2C ) configuration as early infection associated to asymptomatic disease, the ( Figure 2E ) configuration as late infection associated to symptomatic disease, and the ( Figure 2F ) configuration as symbiotic infection, consequent to any approach derived from the application of external driving forces at any time able to maintain configuration ( Figure 2C ) without crashing the system. The goal for any curative approach should be to recover the (Figure 2A) configuration when an infection occurs.

Second, we investigated the system dynamics under different initial conditions, exploring the possible role of different initial viral loads (Figure 3) . Assuming different initial viral loads (10-10,000 RNA copies range), we found a threshold (at about 5,000 RNA copies) for triggering the progressive reduction of Q 1 (Figure 3) . Indeed, for low initial viral load (10-1,000 RNA copies), the system perturbation could be absorbed by the configuration itself (Supplementary Figure 1) , without affecting the overtime stock value of Q 1 and Q 2B but maintaining constant Q 3 , Q 4 , and Q 5 .

The behavior of stock values Q 1 , Q 2A , and Q 2B diverged non-linearly at the threshold value, with a progressive decrease, starting at Day 3 from infection (cyan and red lines, respectively). We found a non-linear evolution of the system output (the Q 5 stock) depending on the initial conditions: for Q 3 stock in the range (10-1,000 RNA copies), the Q 5 was linear, whereas for higher initial viral load, the growth of Q 5 was linear in the first day and non-linear in the further timeframe, associated to a progressive, unpredictable reduction of Q 3 -Q 4 stocks, reflecting in biological terms the turnover of viral proteins required for virions budding.

The non-linear behavior of the (+)ssRNA virus-host interaction was due to the control of Q 4 stock on processes of the host cell, forcing host proteins Q 2B to favor the production of virions (to feed the Q 5 stock value) through the increased J 23 , J 25 , and J 27 flows. In biological terms, the FIGURE 2 | Systems configurations based on initial conditions and effects of external driver forces. In the configuration of initial null viral load (A), the value of stocks Q 1 , Q 2A , and Q 2B were constant, and the system behavior was stationary (B), with constant values of all stocks overtime. At time of infection, the network of flows and feedbacks identified a new configuration (C), to generate a not stationary pattern (D), in which stock values change overtime in response to the other elements of the system, which can evolve to maximize the virions' stock (E). Application of external driving forces, acting on processes (identified by red cross on J 5 ), can reduce the flows and address new configurations (F), identifying leverage points that can be explored at different magnitude and timepoints with a computational simulator.

progressive Q 1 reduction reflects the metabolic rewiring of infected host cells (Chen et al., 2016) , with progressive reduction of resources allocated for the maintenance of host processes, requiring a metabolic shift to less efficient but more rapid source of the available energy required for downstream processes (Thaker et al., 2019 ). The described system dynamics was experimentally validated for influenza virus (Mahmoudabadi et al., 2017) . Results confirmed a previous theoretical assumption, showing a Gibbs free energy for virus lower than its host (Popovic and Minceva, 2020a,b) , when virus and host cells are evaluated separately and not as a unique system as herein proposed. The behavior of Q 1 following different initial values of Q 3 could explain both the variable incubation time in each individual subject and why infections due to (+)ssRNA viruses can occur asymptomatically in most cases. The resilience of virus-host system for a specific range of Q 3 amount could, in turn, depend on intrinsic and extrinsic factors. Indeed, in response to manipulation of stocks and flows, and based on timeframe of observation, the system could evolve along different, non-linear trajectories, requiring early intervention upon infection to make the system resilient to growth of viral stocks.

Currently, the search for a therapeutic strategy is based on single target-related parameters, while we propose to identify systemic targets (i.e., polytarget) to improve host response to host-virus interaction. Starting from the pharmacodynamics of compounds currently under investigation for a typical (+)ssRNA virus (Thorlund et al., 2020) , we could reclassify them, based on their systemic mechanisms of action as listed in Table 3 . Their effects may be potentially simulated to establish the singlecell effect, the best time, and/or schedule of administration, as shortcut of in vitro studies, with a detail level established on the basis of the purpose of the study design. To this end, we applied the search of systemic leverage points by simulating the dynamics of multiple scenarios, upon the action of a generic external driving force (D), assuming that the minimization of the value of Q 5 stock over time should limit the propagation of virions outside the cells.

First, we explored the system configurations upon reduction of the outflows from the Q 5 virions stock, via manipulation of J 7 and/or J 50 . However, minimization of J 7 was counter-effective, due to the increase of Q 5 as consequence of the feedback action in the VR process (data not shown). The effects of full (100%) or partial (50%) reduction of J 50 (flow of energy required for virions budding) could prevent the outflow from Q 5 without stopping its growth, diverting resources from Q 1 and Q 2B to Q 5 , so further supporting viral hijacking of cellular metabolism and impairing host cell homeostasis. As shown in Figure 4 , the effect of driving forces acting to modulate J 50 was different based on application time, Day 0 ( Figures 4A-D) versus Day 1 (Figures 4E,F) , and initial viral load, as the early abrogation of J 50 when the amount of Q 3 was 5,000 RNA copies could restore the stationary status ( Figure 4A) , while halving J 50 maintained homeostasis for hostcell stock, but could not prevent the growth of Q 5 (Figure 4B) . At increasing initial viral load, reduction of J 50 applied at Day 0 (Figures 4C,D) or at Day 1 (Figures 4E,F) could not prevent the growth of Q 5 and the progressive decrease of Q 1 and Q 2B. This systemic dynamics can explain the relationship between the time of initiation of neuraminidase inhibitors and their efficacy (Moscona, 2005) , as treatment starting within the first 12 h after Tenson et al., 2003; Stamatiou et al., 2009; Gielen et al., 2010; Petroni et al., 2020 * Close to Js, refers to azithromycin, in the other colum, being the compount impactinng on the identified flows.

the onset of fever shortened the illness by more than 3 days, as compared with treatment starting at 48 h (Aoki et al., 2003) .

Second, we explored the systemic response to full (100%) or partial (50%) reduction of either J 3 (flow of energy required for RNA replication), or J 4 (flow of energy required for viral RNA translation and viral protein synthesis), or J 5 (flow of energy required for virions assembly), which are involved in virions assembly and are typically dependent on intrinsic viral biological properties. Full reduction of J 5 at Day 0 could recover the systems dynamics in a stationary status, in all scenarios of initial viral load tested (5,000 RNA copies, Figure 5A ; 10,000 RNA copies, Figure 5C ), with constant values for all stocks. A partial reduction of J 5 at Day 0 maintained a stationary status only for lower initial viral load (5,000 RNA copies, Figure 5B ) associated to a reduced amount of stocks Q 4 and Q 5 (Figure 5D) . The abrogation of J 5 at Day 0 was still effective to preserve the stationary status ( Figure 5E ), but halving J 5 at Day 1 could reduce but not prevent the growth of Q 5 (Figure 5F ). We also simulated the effect of applying the same external inputs at different times: after 1 (Supplementary Figure 2) , 3 (Supplementary Figure 3) , FIGURE 4 | System dynamics of (+)ssRNA virus-host interaction in response to external driving forces applied to reduce virions outflow. Changes over time of the values of each stock of the system diagrammed in Figure 1 (for the color code, see bottom), expressed in ATP-eq, in response to reduction of J 50 (flow of energy required for virions budding). Several scenarios are shown: initial viral load 5k and application of full (100%, A) or partial (50%, B) J 50 reduction at Day 0; initial viral load 10k and application of full (100%, C) or partial (50%, D) J 50 reduction at Day 0; initial viral load 10k and application of full (100%, E) or partial (50%, F) J 50 reduction at Day 1.

or 5 days (Supplementary Figure 4) from the initial infection, confirming the role of early application of external driving forces to recover the stationary status of the system.

The importance of full abrogation of J 5 flow for a (+)ssRNA virus-host interaction has been indirectly confirmed by the data recently published by Gordon et al., who cloned, tagged , and expressed 26 of 29 severe acute respiratory syndrome (SARS)-CoV-2 proteins individually in HEK293T cells and used mass spectrometry to measure protein-protein interactions (Gordon et al., 2020b) , to identify 69 existing drugs, known to target host FIGURE 5 | System dynamics of (+)ssRNA virus-host interaction in response to external driving forces applied to reduce virions assembly. Changes over time of the values of each stock of the system diagrammed in Figure 1 (for the color code, see bottom), expressed in ATP-eq, in response to reduction of J 5 (flow of energy required for virions assembly). Several scenarios are shown: initial viral load 5k and application of full (100%, A) or partial (50%, B) J 5 reduction at Day 0; initial viral load 10k and application of full (100%, C) or partial (50%, D) J 5 reduction at Day 0; initial viral load 10k and application of full (100%, E) or partial (50%, F) J 5 reduction at Day 1.

proteins or associated pathways, which interact with SARS-CoV-2, addressing the importance to target the host-virus interaction at the level of RNA translation.

System Dynamics of (+)ssRNA Virus-Host Interaction in Response to External Driving Forces Applied to Reduce Viral Protein Synthesis Full (100%) or partial (50%) reduction of J 4 (flow of energy required for viral RNA translation and viral protein synthesis) at Day 0 did not affect the dynamics of the system (Figure 6) . In particular, for lower initial viral load (5,000 RNA copies, Figures 6A,B) , the Q 5 stock was always lower than Q 2A , thus not affecting Q 1 , which remained constant. However, when Q 5 stock was greater than Q 2A , Q 1 started to decrease, again suggesting that the effect of viral infection on host metabolism is associated to threshold values specific for each system, and not for each cell type. The observation that-based on initial viral loadthe higher value of Q 5 is lower than the stable quantity of Q 2A and Q 2B for 5k ATP-eq and higher than Q 2A and Q 2B for 10k ATP-eq can explain the contribution of initial viral load FIGURE 6 | System dynamics of (+)ssRNA virus-host interaction in response to external driving forces applied to reduce viral protein synthesis. Changes over time of the values of each stock of the system diagrammed in Figure 1 (for the color code, see bottom), expressed in ATP-eq, in response to reduction of J 4 (flow of energy required for viral RNA translation and viral protein synthesis). Several scenarios are shown: initial viral load 5k and application of full (100%, A) or partial (50%, B) J 4 reduction at Day 0; initial viral load 10k and application of full (100%, C) or partial (50%, D) J 4 reduction at Day 0; initial viral load 10k and application of full (100%, E) or partial (50%, F) J 4 reduction at Day 1. and configuration of host-cell stocks in the viral fitness, which depends on host cell cycle stage, addressed as initial value of Q 2A and K 2A (Wargo and Kurath, 2012) .

System Dynamics of (+)ssRNA Virus-Host Interaction in Response to External Driving Forces Applied to Reduce Viral RNA Replication Full (100%) or partial (50%) reduction of J 3 (flow of energy required for RNA replication) at Day 0 did not affect the dynamics of the system (Figure 7) , but-differently from the previous scenario-Q 1 remained constant even if Q 5 was lower than the stable quantity of Q 2A and Q 2B , even if administered later after infection (on Day 3 or 5, as shown in Supplementary  Figures 2-4) . This suggests that the host cell can preserve its homeostasis upon early exposure to inhibitors of viral replication. Interestingly, the pattern in response to single external driving forces was maintained over time, with the values of each stock just shifted, based on the time of application.

Several strategies can be used to reduce selectively virus RNA, including nucleoside analogs, which are metabolized FIGURE 7 | System dynamics of (+)ssRNA virus-host interaction in response to external driving forces applied to reduce viral RNA replication. Changes over time of the values of each stock of the system diagrammed in Figure 1 (for the color code, see bottom), expressed in ATP-eq, in response to reduction of J 3 (flow of energy required for RNA replication). Several scenarios are shown: initial viral load 5k and application of full (100%, A) or partial (50%, B) J 3 reduction at Day 0; initial viral load 10k and application of full (100%, C) or partial (50%, D) J 3 reduction at Day 0; initial viral load 10k and application of full (100%, E) or partial (50%, F) J 3 reduction at Day 1.

intracellularly into their active ribonucleoside 5 -triphosphate forms and incorporated into the nascent viral RNA by errorprone viral RNA-dependent RNA polymerase (RdRps), to disrupt RNA synthesis directly via chain termination, or accumulation of deleterious mutations in the viral genome. The response to nucleoside analogs could be better tested by our computational model, adding the mutations rate and the DNA/RNA metabolism of host cell, not included in the diagram for lack of experimental data about the turnover of RNA stock, specific for each host cell type of interest. 

Third, we next explored the application of multiple external forcing factors, addressing the requirement of a polytarget approach (Bizzarri et al., 2020) . The positive effects on Q 1 , Q 3 , and Q 5 arising from targeting J 5 was mitigated by the combination with reduction of J 50 or J 21 , since the combination targeting J 5 and J 21 (indicated by the blue line in Figure 8 ) applied at Day 0 could preserve the Q 1 amount better than no treatment ( Figure 8A ) associated to slower increase of the growth of Q 3 ( Figure 8B ) and Q 5 (Figure 8C) , leading us to explore further combinations. Application of external driving forces acting on J 50 and J 21 (yellow line in Figure 8 ) increased the Q 5 stock value instead of the expected decrease ( Figure 8C) ; thus, action on a single flow was more efficient than on two of them, given the emergence of compensatory feedbacks and flows. These observations are supported by the controversial findings about the efficacy of macrolides, chloroquine, and their derivatives in the recent COVID-19 pandemic, as their systemic effects include reduction of J 5 and J 21 (Table 3) , and could be affected by the viral load and time of application, with weak changes in virions' spread if applied later in the clinical course of disease as shown by preliminary results of randomized trials (Derwand et al., 2020; Kashour et al., 2020; Magagnoli et al., 2020; Fiolet et al., 2021) .

When the combination targeting J 5 and J 50 was applied later (at Day 5 from infection, orange line in Figure 8 ), Q 1 was higher than untreated ( Figure 8D) , associated to slower growth of Q 3 ( Figure 8E ) and Q 5 (Figure 8F ), but not leading to restore of stationarity. This behavior is a typical systemic feature, where an intervention on a specific local process may lead to counterintuitive rearrangements in the system dynamics. Thus, once the main flow to reduce was found to identify the n-ple stock associated with the desired system output, molecular insights should suggest the biological process to assess a tailored treatment.

In this work, we approached the host-virus interaction dynamics as a systemic problem, and for the first time in the field, we used combined ST tools as a conceptual framework to build up a systemic description of the viral action and host response, critically depending on the existing metabolic environment.

The basic ST idea is to integrate the traditional bottom-up approach-which describes "local" behaviors through causeeffect chains and functional units-with a top-down approach, which points out at the global behavior of the system in terms of its operational configurations, emerging from the feedbackdriven response to different driving forces, like those represented, for example, by the chemistry of new drugs (Odum, 1996) . The utility of computational simulations stems from their capacity to identify structural side effects, non-linearities, and time delays, which are left unexplained by other approaches. Systems biology already recognized the relevance of complexity in the study of microbiological systems (Kaneko, 2006; Loscalzo and Barabasi, 2011) , but although successfully applied in several fields ranging from hard sciences to ecology and economics (Brown, 2004; Brown and Ulgiati, 2004) , the potential of ST in the study of biological systems is still underexploited. Thanks to its abstract nature, stock-flow description can be used in a wide range of different fields, realizing the conceptual bridge that connects the language of biological systems to that of ecology.

In targeting the virus-host interaction, there is an emerging need of tools that could early identify those compounds, not primarily designed for their antiviral action, identifiable by in silico approaches (Gordon et al., 2020b) , which alone or in combination can provide clinical efficacy (Mina et al., 2016; Cheng et al., 2018; Bogdanow et al., 2019; Panja et al., 2019; Gordon et al., 2020b) . There are known advantages of in silico modeling of the action of therapeutic agents on known diseases through agent-based modeling (Mao et al., 2018) . However, the literature evidenced some intrinsic limitations on the choice of parameters, such as the size of investigated populations (Mina et al., 2016) , while major problems are related to model validation (Mina et al., 2016; Donkin et al., 2017) , also requiring to supplement the models with adequate formal ontologies (Kalfoglou and Schorlemmer, 2003; Gotts et al., 2019) .

The proposed model was developed at the single-cell scale. However, in order to define an overall therapeutic approach, the integration with a multiscale approach would be also desirable. In particular, depending on the availability of appropriate data, a future model could focus on different scales, with a more detailed description of some components at the subcellular level, which were grouped (e.g., short-and long-half-life proteins, lipids, and vesicles trafficking) in the present study. On the other hand, the interaction between different cell populations in the host could be also developed, to represent the interaction between healthy and infected cells, and the contribution of immune system (Tan et al., 2007; Katze et al., 2008; Wu et al., 2020) or the repertoire of receptors on the surface of the host cell (Gordon et al., 2020a; Hoffmann et al., 2020; Kailas et al., 2020) to surveil and limit the size of Q 3 stock at the single-cell level. Other natural system constraints could be also included, like some physiological parameters (e.g., temperature and metabolic rate), whose impact on the human body energy dynamics is already understood. The use of a multiscale hierarchical perspective is in principle already possible, as discussed in previous works adopting the same sort of system representation (Ulgiati and Brown, 2009; Ulgiati et al., 2011) .

This work highlights the advantages of applying an ST-based approach to the study of virus-host interaction, being reflected in the possibility of extracting systemic dynamic features that would be otherwise counterintuitive. While a traditional single-target approach would address strategies targeting the viral RNA (Q 3 ) or the replication process (J 3 ), our results suggest that the virus growth is more vulnerable if the process of virion growth before expulsion (process TT, involving flow J 5 ) is targeted.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AR designed the study and collected medical knowledge. FG built up the diagram. MC built up the simulator. All authors performed simulation, analyzed data, prepared the figures, and wrote the manuscript.

This work was supported by SIES, Società Italiana di Ematologia Sperimentale.

The authors thank Luca Naso, Francesco Di Raimondo, and Maria Sanfilippo from SIES, Società Italiana di Ematologia Sperimentale for the support provided.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2020.600254/full#supplementary-material Supplementary Figure 1 | System dynamics of (+)ssRNA virus-host interaction in response to initial viral load. Changes overtime of the values of each stock of the system diagrammed in Figure 1 (for the color code see at the bottom), expressed in ATP-eq: in absence of infection the system status was stationary (A). Upon infection initial values of Q 3 stock identify the response of system, at 50 (B), 100 (C), 1,000 (D), 5,000 (E) and 10,000 (F) RNA copies, expressed in ATP-eq. For Q3 stock in the range (10-1,000 RNA copies) the trajectory of Q 5 evolution was linear, while for higher initial viral load the growth of Q 5 was linear in the first day and non-linear in the further timeframe. 

Mutational and fitness landscapes of an RNA virus revealed through population sequencing

Single molecule analysis of RNA polymerase elongation reveals uniform kinetic behavior

Host Factors in Positive-Strand RNA Virus Genome Replication

Early administration of oral oseltamivir increases the benefits of influenza treatment

Whither systems medicine?

Kinetics of influenza A virus infection in humans

Mengovirus-Induced Rearrangement of the Nuclear Pore Complex: Hijacking Cellular Phosphorylation Machinery

Innate immune evasion strategies of DNA and RNA viruses

Modeling the dynamics of transcriptional gene regulatory networks for animal development

Revisiting the Concept of Human Disease

The dynamic proteome of influenza A virus infection identifies M segment splicing as a host range determinant

A quantitative spatial proteomics analysis of proteome turnover in human cells

Effect of High vs Low Doses of Chloroquine Diphosphate as Adjunctive Therapy for Patients Hospitalized With Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Infection: A Randomized Clinical Trial

Common Nodes of Virus-Host Interaction Revealed Through an Integrated Network Analysis

Systems and models: complexity, dynamics, evolution, sustainability. Norderstedt: Books on Demand GmbH

A picture is worth a thousand words: energy systems language and simulation

Energy quality, emergy, and transformity: H.T. Odum's contributions to quantifying and understanding systems

Evolvability of an RNA virus is determined by its mutational neighbourhood

ATP1A1-mediated Src signaling inhibits coronavirus entry into host cells

A hierarchy of ATP-consuming processes in mammalian cells

Experimental Treatment with Favipiravir for COVID-19: An Open

AMPK regulates energy expenditure by modulating NAD+ metabolism and SIRT1 activity

RNA nuclear export is blocked by poliovirus 2A protease and is concomitant with nucleoporin cleavage

Integrated Pan-Cancer Map of EBV-Associated Neoplasms Reveals Functional Host-Virus Interactions

Six Hours after Infection, the Metabolic Changes Induced by WSSV Neutralize the Host's Oxidative Stress Defenses

Network-based approach to prediction and population-based validation of in silico drug repurposing

Human cytomegalovirus UL37 immediate-early regulatory proteins traffic through the secretory apparatus and to mitochondria

Histone deacetylation is involved in the transcriptional repression of hTERT in normal human cells

Coronavirus NSP6 restricts autophagosome expansion

Global Interactomics Uncovers Extensive Organellar Targeting by Zika Virus

The broad-spectrum antiviral ribonucleoside ribavirin is an RNA virus mutagen

Novel insights into human respiratory syncytial virus-host factor interactions through integrated proteomics and transcriptomics analysis

Current estimates for HIV-1 production imply rapid viral clearance in lymphoid tissues

Cytoplasmic viral replication complexes

COVID-19 outpatients: early risk-stratified treatment with zinc plus low-dose hydroxychloroquine and azithromycin: a retrospective case series study

Viral infection switches non-plasmacytoid dendritic cells into high interferon producers

Antiviral Drug Targets of Single-Stranded RNA Viruses Causing Chronic Human Diseases

Turnover of the human proteome: determination of protein intracellular stability by dynamic SILAC

Determination of the interactome of non-structural protein12 from highly pathogenic porcine reproductive and respiratory syndrome virus with host cellular proteins using high throughput proteomics and identification of HSP70 as a cellular factor for virus replication

Replicating complex agent based models, a formidable task

Proteome half-life dynamics in living human cells

The cellular interactome of the coronavirus infectious bronchitis virus nucleocapsid protein and functional implications for virus biology

HLA class II genes determine the natural variance of hepatitis C viral load

Effect of hydroxychloroquine with or without azithromycin on the mortality of coronavirus disease 2019 (COVID-19) patients: a systematic review and meta-analysis

Dengue virus induces and requires glycolysis for optimal replication

Virus-host interactomes and global models of virus-infected cells

Generation and Comprehensive Analysis of Host Cell Interactome of the PA Protein of the Highly Pathogenic H5N1 Avian Influenza Virus in Mammalian Cells

Elucidation of the Ebola virus VP24 cellular interactome and disruption of virus biology through targeted inhibition of host-cell protein function

System-wide Profiling of RNA-Binding Proteins Uncovers Key Regulators of Virus Infection

Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial

Azithromycin induces anti-viral responses in bronchial epithelial cells

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

Agent-based modelling of socio-ecological systems: Models, projects and ontologies

The Philosophy and Epistemology of Simulation: A Review

A Global Interactome Map of the Dengue Virus NS1 Identifies Virus Restriction and Dependency Host Factors

Extending systems thinking in planning and evaluation using group concept mapping and system dynamics to tackle complex problems

Sensitivity analysis for models with multiple behavior modes: a method based on behavior pattern measures

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

Resolution of primary severe acute respiratory syndrome-associated coronavirus infection requires Stat1

Identification and characterization of homeobox transcription factor genes in Strongylocentrotus purpuratus, and their expression in embryonic development

Interactome analysis of the lymphocytic choriomeningitis virus nucleoprotein in infected cells reveals ATPase Na+/K+ transporting subunit Alpha 1 and prohibitin as host-cell factors involved in the life cycle of mammarenaviruses

Rapamycin suppresses 5'TOP mRNA translation through inhibition of p70s6k

Coronaviruses and the human airway: a universal system for virus-host interaction studies

An interactome map of the nucleocapsid protein from a highly pathogenic North American porcine reproductive and respiratory syndrome virus strain generated using SILACbased quantitative proteomics

Homology Modeling and Docking Studies of TMPRSS2 with Experimentally Known Inhibitors Camostat Mesylate, Nafamostat and Bromhexine Hydrochloride to Control SARS-Coronavirus-2

Ontology mapping: the state of the art

Life: An Introduction to Complex Systems Biology

Efficacy of chloroquine or hydroxychloroquine in COVID-19 patients: a systematic review and meta-analysis

Innate immune modulation by RNA viruses: emerging insights from functional genomics

In vitro inhibition of severe acute respiratory syndrome coronavirus by chloroquine

Characterization of host proteins interacting with the lymphocytic choriomeningitis virus L protein

Antiviral potential of ERK/MAPK and PI3K/AKT/mTOR signaling modulation for Middle East respiratory syndrome coronavirus infection as identified by temporal kinome analysis

A Map of the Arenavirus Nucleoprotein-Host Protein Interactome Reveals that Junin Virus Selectively Impairs the Antiviral Activity of Double-Stranded RNA-Activated Protein Kinase (PKR)

Stochasticity and traffic jams in the transcription of ribosomal RNA: Intriguing role of termination and antitermination

Aberrant innate immune response in lethal infection of macaques with the 1918 influenza virus

Global analysis of host-pathogen interactions that regulate early-stage HIV-1 replication

gative bacteria: the Sec and Tat dependent protein transport pathways

Alteration of protein levels during influenza virus H1N1 infection in host cells: a proteomic survey of host and virus reveals differential dynamics

Interactome Analysis of NS1 Protein Encoded by Influenza A H7N9 Virus Reveals an Inhibitory Role of NS1 in Host mRNA Maturation

Interactome Analysis of the NS1 Protein Encoded by Influenza A H1N1 Virus Reveals a Positive Regulatory Role of Host Protein PRP19 in Viral Replication

Genomic Analysis of the Emergence, Evolution, and Spread of Human Respiratory RNA Viruses

The role of mutational robustness in RNA virus evolution

The Dengue Virus Replication Complex: From RNA Replication to Protein-Protein Interactions to Evasion of Innate Immunity

The incidence of thromboembolism for lenalidomide versus thalidomide in older patients with newly diagnosed multiple myeloma

Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro

Systems biology and the future of medicine

Outcomes of Hydroxychloroquine Usage in United States Veterans Hospitalized with COVID-19

Energetic cost of building a virus

An agent-based model for drug-radiation interactions in the tumour microenvironment: Hypoxia-activated prodrug SN30000 in multicellular tumour spheroids

Exploring the Human-Nipah Virus Protein-Protein Interactome

Thinking in systems : a primer

Cardif is an adaptor protein in the RIG-I antiviral pathway and is targeted by hepatitis C virus

Entrainment and Control of Bacterial Populations: An in Silico Study over a Spatially Extended Agent Based Model

Mathematical model of influenza A virus production in large-scale microcarrier culture

Neuraminidase inhibitors for influenza

Using SILAC and quantitative proteomics to investigate the interactions between viral and host proteomes

Epigenetic Control of rDNA Loci in Response to Intracellular Energy Status

Evaluation of influenza A/Hong Kong/123/77 (H1N1) ts-1A2 and cold-adapted recombinant viruses in seronegative adult volunteers

Histone Deacetylase 2 Is a Component of Influenza A Virus-Induced Host Antiviral Response

Host cell interactome of HIV-1 Rev includes RNA helicases involved in multiple facets of virus production

Viral Innate Immune Evasion and the Pathogenesis of Emerging RNA Virus Infections

Enhancing implementation science by applying best principles of systems science

Environmental accounting : EMERGY and environmental decision making

Explanations of ecological relationships with energy systems concepts

Modeling for all scales : an introduction to system simulation

Modification of drug-binding proteins associated with the efflux pump in MDR-MTB in course of evolution: an unraveled clue based on in silico approach

Molecular Mechanisms of Action of Ribavirin

Human interactome of the influenza B virus NS1 protein

The Synthetic Antiviral Drug Arbidol Inhibits Globally Prevalent Pathogenic Viruses

A complete set of nascent transcription rates for yeast genes

Clarithromycin inhibits autophagy in colorectal cancer by regulating the hERG1 potassium channel interaction with PI3K

The SARS-coronavirus-host interactome: identification of cyclophilins as target for pan-coronavirus inhibitors

Histone deacetylase is a direct target of valproic acid, a potent anticonvulsant, mood stabilizer, and teratogen

Thermodynamic insight into viral infections 2: empirical formulas, molecular compositions and thermodynamic properties of SARS, MERS and SARS-CoV-2 (COVID-19) viruses. Heliyon 6:e04943

A thermodynamic insight into viral infections: do viruses in a lytic cycle hijack cell metabolism due to their low Gibbs energy? Heliyon

Genetic Variability and Evolution of Hepatitis E Virus

Cooperation between translating ribosomes and RNA polymerase in transcription elongation

On the validation of system dynamics type simulation models

Optimal replication of poliovirus within cells

Casein Kinase 2 Is Linked to Stress Granule Dynamics through Phosphorylation of the Stress Granule Nucleating Protein G3BP1

Species-Specific Host-Virus Interactions: Implications for Viral Host Range and Virulence

Role of the Endoplasmic Reticulum-associated Degradation (ERAD) Pathway in Degradation of Hepatitis C Virus Envelope Proteins and Production of Virus Particles

Histone Deacetylase Inhibitor Suberoylanilide Hydroxamic Acid Suppresses Human Adenovirus Gene Expression and Replication

Innate immunity induced by composition-dependent RIG-I recognition of hepatitis C virus RNA

Genome-wide studies in multiple myeloma identify XPO1/CRM1 as a critical target validated using the selective nuclear export inhibitor KPT-276

Global quantification of mammalian gene expression control

Interaction of the innate immune system with positive-strand RNA virus replication organelles

Interferon bioassay: reduction in yield of myxovirus neuraminidases

Broad-spectrum antiviral GS-5734 inhibits both epidemic and zoonotic coronaviruses

Antiviral activity of arbidol against influenza A virus, respiratory syncytial virus, rhinovirus, coxsackie virus and adenovirus in vitro and in vivo

Favipiravir, an anti-influenza drug against lifethreatening RNA virus infections

CX-4945, an orally bioavailable selective inhibitor of protein kinase CK2, inhibits prosurvival and angiogenic signaling and exhibits antitumor efficacy

Transimulation -protein biosynthesis web service

Mechanical and Systems Biology of Cancer

Statistical characteristics of amino acid covariance as possible descriptors of viral genomic complexity

Azithromycin has an antiproliferative and autophagic effect on airway smooth muscle cells

All models are wrong: reflections on becoming a systems scientist

Systems biology and the host response to viral infection

Computational disease modeling -fact or fiction?

The mechanism of action of macrolides, lincosamides and streptogramin B reveals the nascent peptide exit path in the ribosome

Viral hijacking of cellular metabolism

A real-time dashboard of clinical trials for COVID-19

Emergy and ecosystem complexity

Emergy-based complexity measures in natural and social systems

The ER quality control and ER associated degradation machineries are vital for viral pathogenesis

Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population

Chloroquine is a potent inhibitor of SARS coronavirus infection and spread

Rhinovirus 3C protease facilitates specific nucleoporin cleavage and mislocalisation of nuclear proteins in infected host cells

Adjuvant treatment with a mammalian target of rapamycin inhibitor, sirolimus, and steroids improves outcomes in patients with severe H1N1 pneumonia and acute respiratory failure

Comparative influenza protein interactomes identify the role of plakophilin 2 in virus restriction

Host cell interactome of PA protein of H5N1 influenza A virus in chicken cells

Identification of host cellular proteins that interact with the M protein of a highly pathogenic porcine reproductive and respiratory syndrome virus vaccine strain

Viral fitness: definitions, measurement, and current insights

Influenza virus-host interactome screen as a platform for antiviral drug development

Protein turnover

Differential Disruption of Nucleocytoplasmic Trafficking Pathways by Rhinovirus 2A Proteases

Viral dynamics in human immunodeficiency virus type 1 infection

Kinetics of degradation of 'short-' and 'long-lived' proteins in cultured mammalian cells

The mechanism of action of FK-506 and cyclosporin A

Major antigenic site B of human influenza H3N2 viruses has an evolving local fitness landscape

The interactome of the human respiratory syncytial virus NS1 protein highlights multiple effects on host cell biology

Development of an urban FEW nexus online analyzer to support urban circular economy strategy planning

Middle East respiratory syndrome coronavirus ORF4b protein inhibits type I interferon production through both cytoplasmic and nuclear targets

Interactome analysis of the influenza A virus transcription/replication machinery identifies protein phosphatase 6 as a cellular factor required for efficient virus replication

The domain landscape of virus-host interactomes

A Proteomics Survey of Junin Virus Interactions with Human Proteins Reveals Host Factors Required for Arenavirus Replication

Coronaviruses -drug discovery and therapeutic options