Addressing Non-linear System Dynamics of Single-Strand RNA Virus–Host Interaction

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.

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.

Development of a Stock-Flow Diagram
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.

Development of the Virus-Host Interaction Systemic Simulator
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.

Stock-Flow Diagram of (+)ssRNA Virus-Host Interaction
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).

System Dynamics of (+)ssRNA Virus-Host Interaction
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.

System Dynamics of (+)ssRNA Virus-Host Interaction in Response to External Driving Forces Applied to Reduce Virions Outflow
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).

System Dynamics of (+)ssRNA Virus-Host Interaction in Response to External Driving Forces Applied to Reduce Virions Assembly
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.

System Configurations Dynamics in Response to Multiple External Driving Forces Applied at Different Timepoints From Infection
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.

DISCUSSION
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).

CONCLUSION
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.

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

AUTHOR CONTRIBUTIONS
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.

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

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

SUPPLEMENTARY MATERIAL
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.