Executable network of SARS-CoV-2-host interaction predicts drug combination treatments

The COVID-19 pandemic has pushed healthcare systems globally to a breaking point. The urgent need for effective and affordable COVID-19 treatments calls for repurposing combinations of approved drugs. The challenge is to identify which combinations are likely to be most effective and at what stages of the disease. Here, we present the first disease-stage executable signalling network model of SARS-CoV-2-host interactions used to predict effective repurposed drug combinations for treating early- and late stage severe disease. Using our executable model, we performed in silico screening of 9870 pairs of 140 potential targets and have identified nine new drug combinations. Camostat and Apilimod were predicted to be the most promising combination in effectively supressing viral replication in the early stages of severe disease and were validated experimentally in human Caco-2 cells. Our study further demonstrates the power of executable mechanistic modelling to enable rapid pre-clinical evaluation of combination therapies tailored to disease progression. It also presents a novel resource and expandable model system that can respond to further needs in the pandemic.

INTRODUCTION COVID-19 is a complex disease in both dynamics and severity [1][2][3][4] . While most infections only result in mild symptomatic presentation, a proportion suffers more severe disease, for example, the CDC estimate 4.9% of SARS-CoV-2 infections in the USA resulted in hospitalisation in the period up to March 2021 5,6 . Mild cases see an innate immune and interferon (IFN-I/III) response to SARS-CoV-2 infection, as typically observed in other viral infections, likely supporting clearance of the virus and the onset of adaptive immunity 1,7,8 . By contrast, a delayed or absent IFN-I/III response to infection may contribute to the presentation of severe disease 7,[9][10][11][12][13][14][15] . This is especially of concern as emerging variants have adapted to allow potent host immune antagonism 16 , allowing the virus to replicate with little opposition in the early stages of infection 1,10 . This can then be followed by a maladaptively strong and persistent inflammatory response only after the majority of viral replication has occurred 9,17 , leading to life-threatening conditions such as acute respiratory distress syndrome (ARDS) 18 .
There is an urgent need for new effective drugs to manage COVID-19 infection to lower morbidity, mortality and reduce the strain on healthcare systems 19,20 . While several vaccines have been approved [21][22][23] , many millions more patients will need treatment during the years it will take to deliver vaccination worldwide 24,25 . The path for de novo drug discovery is long and complex; the best hope for rapid development of new therapies is, therefore, to combine readily available drugs 26 , with a focus on affordable and well-tested treatments, such as Dexamethasone, which will be necessary to treat COVID-19 in low-income countries that are struggling to obtain sufficient doses of vaccines. Given the large range of potentially suitable compounds, several challenges arise: which drugs are most effective and at what stage of the infection? Are there combinations of drugs that allow effective treatment at lower doses and reduced toxicity?
Driven by the different characteristics of early and late stages of severe COVID-19, and the need to find therapies appropriate to each stage of the disease without interfering with the effective immune response in mild cases, we have developed a computational model that can reproduce these different phases of COVID-19. This approach builds upon our previous success at finding novel drug targets and effective drug combinations for cancer therapy [27][28][29][30] . We created a detailed network map of the interaction between SARS-CoV-2 and lung epithelial cells, as pulmonary involvement, in particular, is a hallmark of severe disease. We use this model to screen thousands of drug combinations to find treatments that block either key viral-host interactions important to viral replication, or the pathogenic dysregulation of the immune response. We focus especially on host-directed therapies, which may be more robust to future variants of the virus 31 . We identify several combinations of repurposed drugs that are predicted to act together to target viral replication in the early stages of the disease or inflammation in the late stage. Moreover, we propose that the model and screen results represent a powerful resource for the community to generate hypotheses about the potential effects of targeting host and viral targets in combination.

RESULTS
Executable network map of the interaction between SARS-CoV-2 and lung epithelium We model the interaction map between SARS-CoV-2 and lung epithelial cells as an executable Qualitative Network 32 model (Fig. 1) using the BioModelAnalyzer (BMA) tool peaked (Fig. 2b, c and Supplementary Table 3). Finally, our model also represents mild disease and uninfected individuals, in order to investigate whether treatments are predicted to have adverse effects beyond known side effects. For example, treatments that reduce maladaptive inflammation in the late severe case may inhibit beneficial inflammatory response if given improperly at the early stages of infection or in mild disease.

Network model validation by comparison to known effects of SARS-CoV-2 infection and treatment
Having generated a network model of viral infection, we first verified that it reproduced key results from experiments that have defined SARS-CoV-2 interaction with lung epithelial cells, as well as known clinical results establishing the utility of different treatments at different stages of the disease (Supplementary Table 4). In particular, we tested the model responses to monotherapies that have been used as COVID-19 treatments, such as the antivirals Remdesivir 35 and Lopinavir 36 ; and the immune modulators Tocilizumab 37 , Dexamethasone 38 , Interferon α-2a (Roferon-A) 39 and Ruxolitinib 40 . Drugs were modelled by setting the affected node in the signalling network to either zero, to represent an inhibitory drug, or to the maximum value, to represent an agonist. The effect of the treatment was then assessed by comparing the resulting level of the nodes representing cellular processes such as viral replication and inflammation to known experimental results (Supplemental Table 4). We then predicted the response of all phenotypes to these drugs (Fig. 3).
Inhibiting protein translation is predicted to suppress viral replication in the early stages of severe COVID-19 Our approach to find new and improved treatments is to test the response of our computational model to an exhaustive search across combinations of many different drugs, either clinically approved or in phase II/III trials, to determine optimal treatments for each stage of development of the disease (Table 1). We , such as a potential treatment, are used as a testing dataset to validate the model. The readouts of the model are compared to the testing dataset and the model is refined iteratively until it reproduces all the experimentally observed behaviours. We then screen the effects of potential drug treatments, either singly or in combination, on the model, to find the best-predicted therapies for early and late stage severe COVID-19. Fig. 2 Computational modelling of host-virus interaction in COVID-19. a The SARS-CoV-2 infection signalling network as modelled in the BioModelAnalyzer tool. Viral proteins depicted in pink, host proteins depicted in green, cellular processes depicted in orange, RNAs depicted in purple. Activating interactions are represented by arrows (→), inactivating by bars (⊣). b Description of the model constraints and outputs for the three states of the disease and the control healthy state. c Trajectories of IFN -I/III response and viral load in mild and severe COVID-19 (as described in Park and Iwasaki, Cell Host and Microbe, 2020) were used to inform the model. implemented the screening by fixing the activity of nodes in the network to a constant value to represent the effect of mono-or combination treatments over all drug combinations and disease stages (Fig. 2b and Supplementary Tables 3, 5). This in silico screen predicted the most effective drug treatments against each stage of COVID-19 and checked for additional adverse effects beyond known drug side effects in mild cases of the disease compared to the severe cases.
At the presentation stage of severe COVID-19, characterised by a weak IFN-I/III response (Fig. 2b and Supplementary Table 3), treatments that are able to arrest viral replication before it can trigger an excessive inflammatory response are needed. As this stage occurs mainly prior to hospitalisation, these treatments must be taken as anticipatory therapy; at or close to the time of suspected exposure or infection. It is therefore critical that they do not adversely affect mild COVID-19 patients, as they will likely be administered before it is possible to stratify patients by prognosis. Our model predicts that, for example, consistent with known treatments in clinical use and trials, Remdesivir 35 and recombinant Interferon α-2a (Roferon-A) 39 are effective at reducing viral replication at the early stage of severe disease (Fig. 4a). In addition, several host translation inhibitors drawn from anti-cancer therapies were also shown to be effective by our model, including Homoharringtonine or Rapamycin (Fig. 4a). Homoharringtonine has previously been identified in a profile of potential drugs targeting SARS-CoV-2 41 , believed to act through inhibition of the eukaryotic small ribosomal subunit (60S) 42 , suppressing viral replication as the virus depends on host translation machinery for the production of viral proteins 43 . Homoharringtonine has cleared Phase II trials for chronic myeloid leukaemia 42 and so is a good candidate to be fast-tracked for use. Similarly, Rapamycin and other mTOR inhibitors, which block translation, are in trials as post-exposure prophylaxis for COVID-19, to prevent the early stage of severe disease from progressing 44,45 . Our model predicts that Rapamycin is only effective against cap-dependent rather than cap-independent protein synthesis. However, as SARS-CoV-2 is dependent upon hijacking the cap-dependent form of protein synthesis 46 this allows equal effectiveness in suppressing viral replication as with other translation inhibitors, with less effect on host protein synthesis and possibly fewer side effects (Fig. 5i). Other translation inhibitors drawn from anti-cancer therapies further demonstrate the potential of host-directed treatments 47,48 .

Combination of Camostat and Apilimod suppress viral replication in early severe COVID-19
Our analysis further reveals options to combine different treatments, focusing on cases where there is an improvement in other pathological processes, in addition to decreasing viral replication. Viral entry inhibitors Camostat (TMPRSS2 inhibitor) and Apilimod (PIKfyve inhibitor) were shown to reduce viral entry but not viral replication in the model, whereas in combination they showed additional ability to prevent viral replication (Fig. 4a, b). This combination also showed improved effects on other key aspects of pathology compared to monotherapy (the effect of mono-vs combination therapy is shown in Fig. 5a for all simulated symptoms), recapitulating the effects of Apilimod seen in TMPRSS2-negative cell lines 49 . Miltefosine (AKT inhibitor) was shown in the model to be ineffective against viral replication alone but exhibits several promising combinations. A combination of Miltefosine with Ulixertinib (MAPK inhibitor) predicted effects across the redundant pathways controlling translation, leading to similar efficacy as Rapamycin (Fig. 4a, b). However, increased inflammation is also predicted as a side effect (Fig. 5b). Conversely, a combination of Miltefosine with the approved anti-inflammatory Dexamethasone was predicted by the model to target both viral replication and inflammation (Figs. 4b, 5c), if the ERK-inhibitory effects of Dexamethasone are sufficient to replace a dedicated inhibitor like Ulixertinib 50 .
Reduction of ER stress using 4-PBA predicted to reduce excessive inflammation in late stage severe COVID- 19 We next considered the late stage of the severe disease, as defined by an inappropriate inflammatory response that must be curbed ( Fig. 2b and Supplementary Table 3). These treatments must be used with care, lest they reduce the useful antiviral immune response in the early severe or mild cases. As seen in the landmark RECOVERY trial 38 , Dexamethasone alone is a partially effective anti-inflammatory treatment for this stage of COVID-19 infection (Fig. 4c). Our model also predicts that sodium phenylbutyrate (4-PBA) may have some efficacy at this stage. Endoplasmic reticulum (ER) stress is dysregulated by the SARS-CoV spike protein (S) 51,52 . It has been suggested that this pathway plays a pro-inflammatory role, via activation of NFκB, in coronavirus infections in general 53,54 ; preliminary results show similar effects in SARS-CoV-2 55 57,58 . However, this therapy is predicted by the model to slightly worsen viral replication in the mild case by decreasing IFN-II activity, and antiviral immunity ( Supplementary Fig. 2h). While this change is small, it may prevent this therapy from being used prophylactically, but it can be safely applied in the late stage of severe disease, once viral replication has peaked, and viral load is waning. Similarly, the IFN-I blocking agent Anifrolumab decreases inflammation in the late stage of severe disease, but this has a counterproductive effect on viral replication in mild disease ( Supplementary Fig. 2h).
Antiviral Lopinavir predicted to suppress inflammation when used in combination with Ruxolitinib in late stage severe COVID-19 Ruxolitinib is already predicted to be effective as a monotherapy and is currently being trialled for COVID-19 40 . However, Ruxolitinib alone may, through inhibition of JAK-STAT and IFN-I signalling,   Table 6.
impair innate immune-mediated suppression of viral entry and replication in mild cases ( Fig. 3f and Supplementary Fig. 2g, h), and has no benefit in early severe disease ( Fig. 3f and Supplementary Fig. 3d). While our model does not suggest combination treatments that can ameliorate this, meaning that Ruxolitinib is likely too risky to use in mild disease, we predict several ways to enhance the effects of Ruxolitinib appropriate for the late stage of severe disease. These include drugs that are ineffective alone such as SAR113945 (IKK inhibitor), Infliximab (TNF-α inhibitor), Canukinumab (IL-1B inhibitor) and recombinant IL-10 (Figs. 4c, d, 5e-j).
The trade-off of these combinations suggests that they should be used only in the late stage of the disease. For example, a combination of Rapamycin and Ruxolitinib is predicted to be more effective against inflammation than either drug alone. However, together these drugs are also predicted to increase viral entry (Fig. 5i). This may suggest that Rapamycin alone could be a good early monotherapy and, if the disease progresses, could be combined with Ruxolitinib in the late stage when a viral entry has already peaked. However, secondary infection, especially bacterial, could be a concern given the potent immunosuppressive effects of some of these combinations. Lopinavir, often administered with Ritonavir, an agent that increases the half-life of Lopinavir in plasma, acts against the virus directly by blocking viral proteases, but has not demonstrated success in treating hospitalised COVID-19 patients 36,59 . This is consistent with our model's prediction, which suggests that while Lopinavir may be effective at reducing viral replication in the early stages of severe COVID-19 (Fig. 4a), it is not predicted to provide additional benefit in the late stage of the disease ( Supplementary  Fig. 4g, h). However, while Lopinavir inhibits the 3CL protease of SARS-CoV-1 60 , it may only be effective against SARS-CoV-2 at toxic doses 41 or not at all 61 though dedicated 3CLPro inhibitors are in development 62 , including Paxlovid (PF-07321332) 63 . Lopinavir is one of the therapies currently being tested in the FLARE trial 64 , which may shed light on this.
Promisingly, targeting 3CL protease in our network model showed interesting anti-inflammatory effects, especially when used in combination with other drugs. 3CL protease is responsible for the proper processing of many viral proteins, including those involved in RNA-dependent RNA polymerase activity (RDRP). Without RDRP activity, the virus is unable to produce sgRNAs required for translation of structural and accessory proteins, including p3a and S, which have been shown to have a proinflammatory effect 55,56 . Direct targets of 3CLpro activity like viral nsp10 have also been shown to modulate the expression of proinflammatory cytokines 65 . This means that 3CL protease inhibition could potentially increase the effectiveness of anti-inflammatories such as Anifrolumab (IFNAR inhibitor) or Ruxolitinib (JAK1/2 inhibitor) (Figs. 4c, d, 5j, k). As such, anti-3CLpro drugs, such as Lopinavir or Paxlovid, may provide benefit even after viral replication has peaked by also playing a role in combination with other drugs to control dysregulated inflammation, and so can be beneficial in both phases of the disease. This contrasts with drugs such as Remdesivir, which blocks viral replication effectively by blocking RDRP but leaves other pro-inflammatory proteins intact. Consequently, our model predicts that Remdesivir will only be effective when applied early in the course of the disease. This is in line with clinical trials showing only modest effects when applied in patients who are already hospitalised 35 .
Characteristic differences in cytokine level distinguish mild, early severe, and late severe COVID-19 As our model predictions show that the effect of drugs is dependent on the stage of the disease at which they are administered, we further searched using our network model for potential prognostic biomarkers for the different stages by comparing the steady-state levels of each node in mild, early, and late severe COVID-19. We observed that a characteristic signature of the early disease is a lower activity of cytokines, such as TNF-α, IL-6 and IL-10 and Interferons type I-III, compared to the mild case ( Supplementary Fig. 5). These will subsequently rise in the late stage of the severe case, with some such as IL-6 exceeding the level seen in the mild case. This is in line with the evidence that severe disease arises, in part, due to an insufficient initial innate immune response, followed by a maladaptively strong and persistent inflammatory response and the prior identification of TNF-α, IL-6, IL-8, IL-10 and CXCL10 as prognostic markers for COVID-19 disease severity in hospitalised patients 33,[66][67][68][69] . This suggests that measuring cytokine levels could help distinguish mild from severe cases in the early stages of the disease, but we also note the risk of an incorrect prognosis due to similarities between late severe and mild cases.

Camostat and Apilimod combined have significantly greater effect suppressing viral entry and replication
To validate our predictions, we selected the most promising combination treatments affecting viral replication identified from modelling early severe disease. Our model predicts that Camostat and Apilimod would have an enhanced effect on viral replication in combination compared to monotherapy (Supplementary Table  7). We tested this in Caco-2 cells and found that, as has been previously reported 70 , Camostat is effective at limiting viral entry leading to a reduction in SARS-CoV-2 nucleocapsid proteinpositive cells in culture, but cannot completely suppress viral entry and replication (Fig. 6a). The addition of Apilimod significantly reduced SARS-CoV-2 infection further, even at the maximum dose of Camostat (Fig. 6b, F = 90.67, p = 2.8 × 10 −10 ), in an additive manner (Supplementary Fig. 6). We further investigated whether the reported inhibitory effects of Dexamethasone on ERK activation 50 were sufficient to synergise with Akt inhibition by Miltefosine (Supplementary Fig. 27) to exert direct antiviral effects, as this may provide extra benefit to Dexamethasone application early in disease progression. This builds on earlier findings that the combination of MAPK and PI3K pathway inhibitors is effective for MERS-CoV 71 , and inhibitors of the PI3K/ Akt/mTOR pathway are effective against SARS-CoV-2 72,73 . However, this combination was ineffective except at doses of Miltefosine that showed toxicity ( Supplementary Fig. 7d), potentially due to limited ERK inhibition by Dexamethasone in Caco-2 cells in vitro. This may be due to the slow nature of the effect of Dexamethasone on ERK, and so beneficial effects may not be evident in short-term cell-culture 50 .

DISCUSSION
Using our computational model, we can predict differential responses to therapy at different disease stages and explore the potential benefits of combination therapies over monotherapies for treating COVID-19 patients to identify the most effective combinations ( Table 2). In this study, we have focussed on screening drugs that are clinically approved or in Phase II/III trials as these have the potential to immediately benefit patients in the current crisis and have had their efficacy validated in clinical practice. We further assume all drugs are maximally effective against all their putative targets to find the broadest array of potential therapies to guide experiments. A specific advantage of our computational approach is that we are able to screen a large number of potential drugs and combinations of drugs in a few hours. We have additionally screened for effective mono-and combination treatments for 36 more drugs of interest that are not yet approved or under trial (Supplementary Table 5 and Supplementary Figs. 8, 9). In addition, we have tested all possible hypothetical interventions against all the viral and host proteins in the network modelscreening a total of over 9,000 combinations. These targets, if they can be made druggable, suggest potential novel viral or host-directed targets ( Supplementary Figs. 10, 11). Together, the network model and the drug screening algorithm provide a valuable resource that can be leveraged by the scientific community to generate hypotheses about the effect of potential therapies on cellular processes and complement existing screens of monotherapies 74 .
There are several computational approaches that have been applied to drug screening that are being applied to COVID-19, each with their own strengths and weaknesses. For example, as the structures of SARS-CoV-2 proteins have been determined, existing databases of small molecules have been leveraged to survey thousands of candidates to bind to these proteins 75 . Such studies helped identify 3CLpro as a common target, as well as the use of drugs such as Lopinavir and Remdesivir 26 . Another approach is the use of protein-protein interaction networks, which are particularly suited to identifying host-directed therapies 76,77 . These networks can be analysed for their structure alone and can cover a broad range of potential targets. These have been used to identify, for example, the potential for cancer drug repurposing 47,48 . However, static-network analyses rely on the assumption that certain features of the network structure can identify the best targets, e.g. proximity to disease-associated nodes, but they cannot predict specific effects. Moreover, they rely on pre-existing network databases or require additional curation 78,79 or machine learning 80 for new network types. Mathematical models such as Ordinary Differential Equations can interrogate the dynamics of the disease, but require precise data to fit model parameters, and so can only handle a smaller set of variables 81,82 .
By contrast, our approach combines scalable, executable modelling with transparent, biologically plausible explanations 30,83 since each of our predictions is derived from biological interactions that can be explained in the context of the overall model, which itself is derived from, and consistent with, experimentally verified observations (Supplementary Table 1). The mechanistic explanations from our model advance the fight against the COVID-19 pandemic by increasing our understanding of why some treatments work and others fail. Our model and screening are readily accessible and are built using the BMA opensource and freely available toolset (https://biomodelanalyzer.org) and can be updated as new SARS-CoV-2 variants emerge that may exploit different mechanisms to enter cells and replicate or suppress host defence mechanisms 16 .
We deliberately chose to focus on the innate, rather than the adaptive immune response, specifically in lung epithelial cells. We further prioritised predicting the best targets for optimised drugs, rather than attempting to model the full pharmacodynamics of specific compounds. This allowed us to survey a broad selection of potential therapies for a critical cell type in the severe form of COVID-19. It also allowed us to focus specifically on two scenarios for the innate immune response: first, a low response seen in many patients 9-14 that we believe characterises the early stage, and second, a persistent and excessive response in the late phase of the disease 17,84 .
We considered this broad scope appropriate at this stage of the pandemic, when there is a need for rapid development of new treatments, but without compromising safety. In the absence of computational screening of the kind we advocate, there has been a focus on only a small number of drugs 85,86 , many screened through expensive Phase III trials 86 with high-profile safety concerns in the case of Hydroxychloroquine 87 , one of the most studied 85 . Our computational approach can help accelerate the process of screening more drugs and rejecting those with safety concerns, rapidly guiding clinicians to the most likely candidates.  R. Howell et al.
Exemplifying this, amongst the most promising targets for early severe COVID-19 predicted by our model were Camostat and Apilimod, primarily targeting TMPRSS2 and PIKfyve respectively. We show that these drugs are significantly more effective together than alone in suppressing viral entry and replication in Caco-2 cells (Fig. 6). Combining these drugs effectively blocks the two key pathways for viral entry, via the cell membrane exploiting TMPRSS2 and via the endosomal route, without which the virus is severely limited in routes into the cell 70 . This novel combination builds on prior work blocking cathepsin directly 70,88 and demonstrates its effectiveness with two phase II tested drugs. This combination also limits the range of host target cells the virus can infect. Camostat alone will only be effective in cells that allow TMPRSS2 entry, while the TMPRSS2-negative cell population also targeted by the virus remains permissive 89 . The addition of Apilimod may prevent this by targeting the endosomal entry pathway. While current variants of SARS-CoV-2 primarily exploit the TMPRSS2 entry route 90 , this may be beneficial if the virus evolves to exploit the endosomal route more aggressively. Hostdirected therapies such as this and others suggested by the model may be more resistant to future mutations in SARS-CoV-2 31 .
Our computational model could be developed further to address patient stratification and correct dosing. Many treatments, including those we suggest, depend upon being administered at a specific stage of the disease, making it vital to accurately stratify patients, and necessarily includes factors outside the scope of our model, such as when it is practical to administer treatment. As an example, Remdesivir and Interferon α-2a (Roferon-A) have similar effects in the model, but Interferon α-2a may be significantly easier to administer as it is inhaled, and so may be more appropriate in the early case, compared to Remdesivir that requires intravenous administration. The likely side effects of treatments must also be considered, though it is hoped that treatment courses will be relatively short; Miltefosine with Ulixertinib is likely to produce diarrhoea and vomiting, which may be a severe burden in an unwell patient and pose aerosol risks to medical staff attending to them. Many of these adverse effects are driven by interaction outside the scope of the current model and this could be a fruitful area of expansion. Our network model helps address stratification by identifying potential biomarkers for different stages and severities of disease. It could also be extended to address dosing for both mono-and combination therapies. Finally, as we understand more about the progression of the disease and the set of traits that determine mild versus severe disease, the model could be expanded to better and, in more detail, explore the stages of severe COVID-19.
In this study, we have demonstrated that computational modelling has the ability to rapidly screen and predict new and improved therapies for previously unknown and life-threatening conditions. In particular, through this approach, we have listed several new combination therapies, based on existing and approved drugs, with the potential to improve outcomes for COVID-19 patients. The flexibility and transparency of mechanistic computational modelling will allow this work to be further developed as the virus evolves and demands further changes in therapeutic practice.

Qualitative networks
We model the viral-host interaction network as a discrete qualitative network 32 . Qualitative networks are an extension of Boolean networks 91 in which each node may take multiple finite values in a fixed range (e.g. 0,1,2) rather than only ON and OFF. We build and analyse this network model using the open-source (MIT License) and freely available BioModelAnalyzer (BMA) tool 92 https://biomodelanalyzer.org. The model is available at https://github.com/JFisherLab/COVID19. The network consists of nodes representing viral proteins and RNA; host genes and proteins; processes such as viral entry; and phenotypes such as inflammation (see Supplementary Notes). The interactions between nodes are represented as edges (see Supplementary Table 1). The level of activity of a node is represented as a positive integer within a fixed range specific to that node. The level of activity changes in response to other nodes, as determined by a mathematical function associated with the node called a target function. The target function takes as input the level of neighbouring nodes that have an incoming edge to the node the target function controls. The default target function is avg(pos)-avg(neg). This takes the mean of the level of activity of all the nodes with a positive edge (represented by an arrow (→)) and subtracts the mean of the level of activity of all the nodes with a negative edge (represented by a flat head (⊣)). More complex functions are needed to describe nodes with behaviour such as only activating when an input rises above a threshold. These target functions, and their rationale, are described in Supplementary Table 2. If node X has a granularity of a-b and is used in the target function of another node X' with range a'-b', then it is scaled to the range of X′ using the equation: For a given set of initial values for all nodes, the network model will update all node values synchronously, and so there will be a single stable attractor the network will tend to. This attractor may be of any number of states, if it is one state we refer to it as a fixed-point attractor, if greater than one state we refer to it as a loop. We can test for whether the network reaches a fixed-point attractor for all possible initial states 93 . If there are different attractors for different initial states, we refer to this as a bifurcation. In the case of a fixed point, we take the level of every node in the attractor as the prediction of biological behaviour, in the case of a loop or a bifurcation we use the mid-point of the upper and lower bound that can be placed upon node behaviour using the algorithm described in Cook et al. 93 . Details of the procedure of network simulation can be found in Schaub et al. 32 and the bounds and reported mid-point are given in Supplementary Data S1.

Model testing
We built the network using a bottom-up iterative approach. We collated data from the literature describing interactions between and within SARS-CoV-2 and the host cell. The experimental evidence for each edge is described in Supplementary Table 1. We then collated a separate set of testing data used to evaluate the model. This was formed from published data on the effect of drugs on SARS-CoV-2 in culture and is summarised in Supplementary Table 4. These data describe experiments that do not define a point-to-point interaction between genes and/or proteins, but rather the overall behaviour of the system that emerges from the sum of the interactions and so provides a testing dataset for the model separate from the data used to build it.
In order to reproduce an experimental condition in the model, we set the target function of the relevant node to a constant value. For example, the presence of SARS-CoV-2 is represented by changing the target function of Virion to 1. We then find the stable states of the model in this case and compare the predicted stable value of nodes for which there are data (e.g. inflammation or level of IL-6) to that observed in the experiment. We iteratively tested the model and based on these tests, expanded and refined the model until it matched all the observed behaviours, as described in Supplementary Table 4.

In silico drug combination screens
In order to find the most effective treatments, we inactivate (set target function to a minimum) or activate (set target function to maximum) all nodes, or all sets of nodes targeted by a corresponding drug (Supplementary Table 5) either singly or in a pair-wise combination, using the BMA Command Line tool BioCheckConsole. We consider all combinations of nodes and all combinations of drugs, where a single drug may target multiple nodes. This assumes that drugs are able to act on all their putative targets and does not account for other pharmacodynamic effects that may reduce their efficacy. For each node in the network, we manually curated from the literature a list of drugs that could target nodes within the network, finding 74 candidate drugs (Supplementary Table 5), of which 38 had passed phase II trials at the time of writing (Table 1). We also found all the reported targets for these drugs using the DrugBank database 94 to ensure all targets included in the model are direct and clinically relevant. We evaluate these perturbations by the stable state of the network for the eight phenotypes of interest, or in the case that the steady-state could not be found, the mid-point of the upper and lower bounds BMA could place on the behaviour of that node using the algorithm described in Cook et al. 93 . We compare different stages of the disease by using different backgrounds; setting certain nodes to a constant level to reproduce the mild, early severe, and late severe forms of the disease (see Supplementary Notes, Supplementary Table 3 and Supplementary Data 1).

Inhibitor treatment and infections
Caco-2 cells were pre-treated with Camostat Mesilate (ApexBio, B2082), Apilimod (Selleckchem, S6414), Miltefosine (Cayman Chemicals, 63280), Dexamethasone (EMD Millipore, 265005) or DMSO (Sigma) vehicle control at 2x the final concentration for 2 h prior to infection in 50 μl culture medium. Cells were infected with 1000 E copies/cell SARS-CoV-2 in 50 μl, bringing the total culture volume to 100 μl and inhibitors to 1x final concentrations as indicated. Cytotoxicity of inhibitors was determined by tetrazolium salt (MTT) assay. About 10% MTT was added to culture media and cells were incubated for 24 h at 37°C. Cells were lysed with 10% SDS, 0.01 M HCl and the formation of purple formazan was measured at 570 nm.

Flow cytometry
Infection levels were measured at 24 h by flow cytometry. Caco-2 cells were trypsinised, stained with fixable Zombie UV Live/Dead dye (BioLegend) and fixed with 4% PFA before intracellular staining for nucleocapsid protein. For intracellular detection of SARS-CoV-2 nucleoprotein, cells were permeabilised for 15 min with Intracellular Staining Perm Wash Buffer (BioLegend). Cells were then incubated with 1 μg/ml CR3009 SARS-CoV-2 cross-reactive antibody (a kind gift from Dr. Laura McCoy) in permeabilisation buffer for 30 min at room temperature, washed once and incubated with secondary Alexa Fluor 488-Donkey-anti-Human IgG (Jackson Labs). All samples were acquired and analysed using a NovoCyte (Agilent) and NovoExpress 1.5.0 software (Agilent). Supplementary Fig. 13 shows the gating strategy applied to a representative sample of Caco-2 cells.

Calculation of drug combination indices
The expression of intracellular SARS-CoV-2 nucleocapsid protein was measured by flow cytometry at 24 h post infection and used as a measure of drug effect on viral entry and replication. The data were normalised to the value of the control group and averaged across nine replicates. For each treatment, the combined effect of drug 1 at concentration a and drug 2 at concentration b was calculated as E a;b ¼ ð100 À N a;b Þ=100, where N a;b is the average percentage of cells expressing nucleocapsid protein at this concentration. The combination index (CI) was then calculated according to the Bliss independence model 95,96 : When CI < 1, this indicates synergy between the drugs at this concentration, while CI % 1 indicates additivity and CI > 1 indicates antagonism.

Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY
Authors can confirm that all relevant data are included in the paper and/or its supplementary information files. The model is available at https://github.com/ JFisherLab/COVID19.