Communication network model of the immune system identifies the impact of interactions with SARS-CoV-2 proteins

Interactions between SARS-CoV-2 and human proteins (SARS-CoV-2 PPIs) cause information transfer through biochemical pathways that contribute to the immunopathology of COVID-19. Here, we present a communication network model of the immune system to compute the information transferred by the viral proteins using the available SARS-CoV-2 PPIs data. The amount of transferred information depends on the reference state of the immune system, or the state without SARS-CoV-2 PPIs, and can quantify how many variables of the immune system are controlled by the viral proteins. The information received by the immune system proteins from the viral proteins is useful to identify the biological processes (BPs) susceptible to dysregulation, and also to estimate the duration of viral PPIs necessary for the dysregulation to occur. We found that computing the drop in information from viral PPIs due to drugs provides a direct measure for the efficacy of therapies.


INTRODUCTION
COVID-19, the disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has multiple routes for harming human biology. One mechanism is through protein-protein interactions (PPIs) between SARS-CoV-2 and human proteins (SARS-CoV-2 PPIs) [1]. The viral PPIs can alter biochemical functions of human proteins, the effect of which can propagate through the cascade of biochemical reactions informing the abundance of other proteins that do not directly interact with the viral proteins. The information transfer caused by SARS-CoV-2 PPIs may elucidate some of the more complex consequences of COVID-19, like long lasting immunological dysfunction [2][3][4], dysregulation of myeloid response [5], lymphopenia [6,7], and more. For a comprehensive assessment of COVID- 19, we need to quantify the information about SARS-CoV-2 PPIs that is communicated to the immune system proteins and the sensitivity of biological processes (BPs) to the amount of transferred information. The fundamental framework for communication is information theory [8], which has been used to quantify the efficacy of biochemical signaling [9][10][11]. Information theory uses entropy to measure the dependence between stochastic variables, and all biochemical entities are stochastic at the cellular level. By bringing the fundamental information-theoretic model of biochemical signaling to the immune system, we can computationally determine the ramifications of COVID-19 and other infections. Computing the information transfer in the immune system due to viral PPIs also provides a framework for in silico evaluation of therapies [12].
Here, we developed a communication network model of the immune system using pathways and interaction data to compute the information transferred by viral PPIs, and to subsequently identify potentially dysregulated BPs and network updates the state of the network as where X (k) is the network state after k steps of information propagation. If the initial state is X (0) , then after k steps of information transfer the state is X (k) = C k G X (0) . The stationary state of the network satisfies X = C G X. There are 89 source nodes in the network (i.e., nodes that only send but do not receive information), the state of these nodes determine the reference state of the immune system communication network in the absence of other information sources, like SARS-CoV-2 PPIs. Current limitations of the communication network model include, incompleteness and errors in the pathways and interactions data, a binary channel to model the reactions and interactions, and unbiased interference of information. Proteomic, metabolomic, and rate constants data for human cells can provide data-driven parameters for channel bit size and interference bias. We expect to progressively integrate such data in the future versions of the model.

INFORMATION TRANSFERRED BY SARS-COV-2 PPIS
To evaluate the information transferred by SARS-CoV-2 PPIs, we assumed that interactions with SARS-CoV-2 proteins diminish the function of human proteins. So, the impact of SARS-CoV-2 PPIs is a reduction in the functional or the effective abundance of the associated human proteins. For each of the proteins that directly interact with SARS-CoV-2 proteins (SARS-CoV-2 PPI node), we constrained the state of that protein node to {0, 1}, or effectively to low abundance. We need experimental studies to assess the exact decline, if at all, in a protein's function due to SARS-CoV-2 PPIs (as available for interactions with the ORF9b protein [18]). However, we can computationally test all possible impact on the protein function, i.e., maximum decline in function (constraint {0, 1}), no effect on the function (no constraint), or maximum increase in function (constraint {1, 0}), or any intermediate constraint between {1, 0} and {0, 1}. We used the constraint {0, 1} to compute the consequence of maximum decline in a protein's function due to SARS-CoV-2 PPIs.
To quantify the information transferred by SARS-CoV-2 PPIs, we first computed the reference stationary state X ref in the absence of the SARS-CoV-2 proteins, and then constrained the SARS-CoV-2 PPI nodes as described above to compute a new stationary state X using Eq. (1) (SI: 3). The information transferred to a single network node v i due to the SARS-CoV-2 PPIs is given by the relative entropy (or the Kullback-Leibler divergence) in bits as [8] The network relative entropy H ref (X) is the sum of the relative entropy of all the unconstrained network nodes. We computed the information transferred due to SARS-CoV-2 PPIs for multiple reference states (Fig. 1b) BPs as significant (Fig. 2), but also shows that more than one SARS-CoV-2 protein can influence the same BP due to information transfer (e.g. immune effector process and cell activation in Fig. 2). The effect of COVID-19 on cytokine production and inflammation is well known [22,23]. But several of the other BPs that are significant due to communication have also been reported. These include, rapid decay in humoral immunity [24,25], increased apoptosis in infected cells [26,27], localization of organelles, especially endoplasmic reticulum, endosome, and lysosome [28,29], modulation of catabolic recycling processes [30], and protein import into the mitochondria [31], to name a few. Additionally, dysregulation of some of the cell activation processes due to COVID-19, release of myeloid dendritic cells [5,32] and lymphocyte production [6,7], have also been reported. Protease inhibitors, like nirmatrelvir/ritonavir (sold as Paxlovid), reduce viral replication and consequently decrease the duration of SARS-CoV-2 PPIs. Therefore, using the Q-N curves we can identify the BPs that can get dysregulated due to any delay in the use of protease inhibitors.

EVALUATION OF CANDIDATE THERAPIES
To use the communication network model for evaluating candidate therapies, we computed the reduction in the information transfer from SARS-CoV-2 PPIs due to drugs. We demonstrated this technique using the maximum entropy reference state (as in Protocol 1 in Fig. 1b) and then applying the constraint {0, 1} on all SARS-CoV-2 PPI nodes together (all proteins in Table SI: 2), which produces a total network relative entropy of 65.4 bits. The Reactome database contains 19 drugs that interact with the immune system, majority of which are immunomodulators [33]. We constrained each of the drugs to the state {1, 0}, one at a time, and computed the subsequent change in network relative entropy, ∆H ref (X) (Fig. 4a). The maximum net reduction in the SARS-CoV-2 PPIs signaling was 2 bits (for ibrutinib and acalabrutinib) and the maximum net increase was 3 bits (for omalizumab). Or, at best some of the drugs can effectively free 2 out of the 65 variables of the immune system that are under the control of SARS-CoV-2 PPIs. 12 out of the 19 drugs produced less than ±1 bit change in the network relative entropy, indicating that majority of these 19 drugs have minimal effect in countering the signaling from SARS-CoV-2 PPIs. However, a drug can reduce the information transfer from SARS-CoV-2 PPIs in some nodes of the immune system (drop), but it can also cause additional information transfer in other nodes of the immune system (gain). So, for a more nuanced evaluation of the efficacy of a drug, we identified the network nodes that have a drop in relative entropy due to a drug and in which nodes there is a gain, to obtain the cumulative drop/gain values (Fig. 4b which inhibit cytokine activity, were considered to be safe [34,35], but their safety profile is being re-evaluated [36]. Omalizumab, an anti-Immunoglobulin E drug, has several side effects, including anaphylaxis [37]. JAK inhibitors have the best drop values, and they have been shown to increase survival during severe COVID-19 in clinical trials, but insignificant for reducing hospitalization [38][39][40]. Bruton's tyrosine kinase (BTK) inhibitors, which impede the B-cell receptor pathway, have the second best drop values, and ibrutinib has been claimed to reduce pulmonary injury due to COVID-19 in a clinical trial [41]. Interleukin 6 (IL-6) inhibitors emerge as comparatively less effective in our analysis, which is in accordance with the results of clinical trials -IL-6 inhibitors are significantly beneficial only when administered with glucocorticoids in clinical trials [42,43]. But the limitations in our model is probably underestimating the efficacy of tocilizumab, which is one of the recommended treatments for critical COVID-19 [44].
As for the other drugs, the relatively low efficacy of Interleukin-1 (IL-1) inhibitors agrees with the clinical data [45,46].
Programmed cell death protein 1 (PDCD1) inhibitors, which are immune checkpoint inhibitors, are hypothesized to treat COVID-19 induced lymphopenia [47], but we did not find any study to support the claim. Omalizumab is still in a clinical trial for COVID-19 treatment [48]. The current NIH treatment guideline for COVID-19 recommends only baricitinib and tofacitinib among the JAK inhibitors and mostly against the BTK, IL-6, and IL-1 inhibitors, except for the cautious use of tocilizumab and no recommendation for or against anakinra [33].
We similarly computed the effectiveness of artificially controlling the expression level of proteins to reduce the information transfer from SARS-CoV-2 PPIs. One at a time, we set the immune system proteins (except the SARS-CoV-2 PPI nodes) at high abundance {1, 0} and computed the subsequent change in the network relative entropy.
We used artificial overexpression of proteins in this study because it is opposite to the effect of SARS-CoV-2 PPIs on protein function in our model. We found that 638 out of the 1122 proteins that we surveyed produced a better ∆H ref (X) than any of the 19 drugs (i.e, more than 2 bits, Fig. 4c). In particular, we found 23 proteins that can reduce the signaling from SARS-CoV-2 PPIs by 28 to 30 bits, and this set contains 3 IFITM (interferon-induced transmembrane) proteins. Interestingly, artificial overexpression of IFITM proteins has been shown to block SARS-CoV-2 infection in vitro, even though endogenous overexpression of IFITM proteins supports the infection [49]. PSMB8 (Proteasome subunit beta type-8 protein, involved in antigen processing) has been reported to cause milder infection under endogenous overexpression using data from lung samples [50]. Several of the other 23 proteins identified in this analysis (e.g. MUC1, MX1, MX2, IFIT1, XAF1, EGR1) have been reported to be overexpressed during COVID-19, depending on the severity of the disease [51][52][53][54][55]. We separated the drop/gain contributions to ∆H ref (X) for each of the proteins (Fig. 4d). The 23 proteins that produced the best ∆H ref (X) remain close to the top hundred even after ranking the proteins by the drop component. So, it is possible to avoid the proteins with high gain values without compromising on the reduction to SARS-CoV-2 PPIs signaling. Since controlled expression of the proteins listed in Fig. 4c can reduce the information transfer considerably better than any of the drugs, we may be able to create more effective treatments for SARS-CoV-2 PPIs compared to the ones in Fig. 4b.

CODE AVAILABILITY
The source code to construct the communication network from the Reactome and the Interactome databases and to perform the subsequent computations is available at github.com/sarkar-s/core. The enrichment analysis was done using the GOBP gene sets in the MSigDB database (SI: 5) [56,57]. The significant BPs are labeled using their adjusted p-value (Q) after multiple hypothesis testing. The first block of Q values (shaded gray) is for the PPI set for each SARS-CoV-2 protein, and the remaining three blocks are for the sum of the PPI and communicated sets. We identified the significant BPs using the hypergeometric test and a false-discovery rate of 1% [58], which is the same value used by the authors of the SARS-CoV-2 protein interaction map [1]. The top 3 significant BPs for each SARS-CoV-2 protein in Protocol 2 were selected for this figure (top BPs), because Protocol 2 produces the highest amount of information transfer in the network. We show the Q values only for the top BPs using the proteins from the PPI set, Protocol 1, and Protocol 3, to compare the effect of information transfer. The complete list of all significant BPs for each protocol is in SI: 5.  To construct the communication network model of the immune system we used the Reactome biological pathways [1] and the Human Reference Interactome (HuRI) [2] databases. The Reactome database contains the list of biochemical reactions within the top level pathway of the Immune System. In this database, the immune system contains 1622 biochemical reactions and 176 regulatory reactions. Each reaction has a set of inputs and outputs, and each input (or output) is a PhysicalEntity in the Reactome database. Subclasses of PhysicalEntity are EntitiesWithAccessionedSequence (for proteins), Complex (for complexes), SimpleEntity (for molecules like ATP and ions like Ca 2+ ), and some other classes with a smaller number of physical entities. We created the immune system communication network G = (V, E), by adding the edges (or connections) between the PhysicalEntities using the input-output relationships from the ReactionLikeEvent class in the Reactome database. The ReactionLikeEvent class contains the subclasses of Reaction, BlackBoxEvent, Polymerisation, and Depolymerisation. Edges from the Reaction subclass are bidirectional, and the edges from the remaining subclasses are unidirectional. Information (or stochasticity) from each input can propagate to the information in all the outputs of an event. Therefore, for a single event, each pair of input and output is represented by an edge (or an information channel). Hence, a reaction like, A + B −−→ C + D + E, contributes 6 edges to the immune system communication network, {A → C, A → D, A → E, B → C, B → D, B → E}. Inhibitory and negative regulation reactions reduce the abundance of the reaction output with increasing abundance of the input (or the regulator), we represent these edges as A − | | B. All other types of reactions, where the input increases the abundance of the output, are represented as A → B. After constructing the network from the Reactome database, we used the HuRI database [2] to identify the pairs of proteins that interact with each other. Each pairwise interaction was added as a bidirectional edge (type →) to the network (currently, the HuRI database does not contain the sign of interactions). We used the NetworkX library in Python to create the network object. Once the network was created we removed all isolates (i.e., nodes that are not connected to any other nodes) from the network. The total number of nodes in the network is |V | = 3011, and the breakdown of the number of nodes by their class is in

SI: 2. INFORMATION TRANSFER THROUGH THE NETWORK
To simulate the transfer of information in the immune system communication network, we modeled each edge of the network as a binary information channel. Reactions where the abundance of the input increases the abundance of the output (including positive regulation reactions) were modeled as v i → v j , with the input and the output state related as where ρ < 0.5 is the channel error that reduces information transfer from the input v i to the output v j . Reactions where the abundance of the input decreases the abundance of the output (inhibitory and negative regulation reactions) were modeled as bit flip information channels, represented as v i − | | v j , with the input and the output state related as The matrices in equations (1) and (2) are transition matrices for the information channel model for individual edges in the network. We built the global transition matrix for the entire immune system network using the method described in [3]. Let the adjacency matrix of the graph be A, which is a |V |-by-|V | matrix, with A[i, j] = 1 if there is an edge where deg − (v j ) is the in-degree of node v j , or the node that is receiving information from the source v i . The weights w ij determine the type of interference of information. For unbiased interference all weights are w ij = 1. For biased interference the weights has to satisfy the conditions 0 < w ij and j w ij = deg − (v j ). We used unbiased interference in our calculations. Proteomic and metabolomic profile data of human cells [4] can be used to create a data-driven biased interference model. We also used ρ = 0 for the results presented in the main text. So, our results capture the information loss due to signal interference rather than communication noise. Non-zero communication error will further decrease the network relative entropy from the SARS-CoV-2 PPIs, communicate to a smaller set of proteins, and thereby affect a smaller set of biological processes. Data on the rate constants of human proteins and transcripts [5][6][7] can be used to estimate the ρ value for each reaction [8].
We represent the global state of the communication network using a column vector X with the dimension 2|V |. As mentioned in the main text, the i th 2D block of X contains the state of the node v i , {P (v i = 1), P (v i = 0)}. The transfer of information through the network is computed using the equation where each multiplication by C G causes information to transfer from a network node to its nearest neighbors. The stationary state of the network is obtained when X (k+1) − X (k) 1 < , where · 1 is the L1-norm and the is an error tolerance, we used = 10 −4 in all of our calculations, which we found sufficient to achieve convergence in the network relative entropy for the 3 protocols. Finally, it is not necessary to model each node as a binary variable and each edge as a 1 bit information channel. We can easily construct the global transition matrix C G using a n > 1 bit channel model for each edge, using the procedure described through equations (1), (2), and (3). In general, the state vector X has a dimension 2 n |V | and the global transition matrix is 2 n |V |-by-2 n |V |, as a function of the bit size for the channel model for each edge.

SI: 3. INFORMATION TRANSFER DUE TO SARS-COV-2 PPIS
To compute the information transfer in the immune system network due to SARS-CoV-2 PPIs, we first identified the immune system proteins that directly interact with each SARS-CoV-2 protein using the protein interaction map published by Gordon et. al. [9]. The set of immune system proteins that directly interact with each of the SARS-CoV-2 proteins are in Table. SI: 2. -2 proteins PPI set of proteins  NSP2  SLC27A2, EIF4E2  NSP7  CYB5R3, RALA, RAB5C, RAB7A, RAB10, RAB14, RHOA, PTGES2, RAB18  NSP8  HECTD1  NSP10  AP2A2  NSP12  RIPK1  NSP13  RDX, ERC1, TBK1  NSP14  GLA, IMPDH2  NSP15  RNF41  Spike  GOLGA7  ORF3a  HMOX1  E  SLC44A2  M  STOM, ANO6  ORF8  ITGB1, PVR, NPC2, GGH, IL17RA, NEU1, ERP44  ORF9b  TOMM70  ORF9c  NLRX1, ECSIT  N  CSNK2B  ORF10 ELOC, ELOB To determine the information transfer due to each SARS-CoV-2 protein, we constrained the corresponding set of immune system proteins from Table SI: 2 to the state {0, 1} and computed the network relative entropy of the stationary state with respect to the reference state for each protocol. We used ρ = 0 for the results shown in the main text. The network relative entropy values for the 3 protocols shown in Fig. 1(b) are in the file SIdata3. We present the network relative entropy values for a set of other ρ values in Fig. SI: 1. A higher communication error reduces the information transfer, which means a smaller number of immune system proteins will receive information about the SARS-CoV-2 PPIs. Network relative entropy always decreases with increasing ρ. But ρ > 0 can also change the ranking of the SARS-CoV-2 proteins based on the relative entropy induced by them.

SI: 5. SIGNIFICANT BIOLOGICAL PROCESSES DUE TO THE COMMUNICATED SET OF PROTEINS
To identify the biological processes that are candidates for dysregulation, we performed Gene Ontology (GO) overrepresentation analysis using the Biological Processes (BP) gene sets [10]. We used the GOBP gene sets in the c5 collection of the MSigDB database [11]. We determined the total interaction set for each SARS-CoV-2 protein, i.e., sum of the PPI set (SI: 2) and one of the communicated sets (SI: 3-SI: 5) depending on the protocol. Then we determined if the total interaction set is overrepresented in any of the GOBP gene sets. The p-value for each GOBP term was computed using the hypergeometric survival function from scipy.stats in Python. A false discovery rate (FDR) of 1% was used for multiple hypothesis testing using the Benjamini-Hochberg correction, which provided the adjusted p-value (Q). This computation was done using the fdrcorrection function in the statsmodels library. After determining the significant GOBP terms (for FDR 1%), we further reduced the set of significant terms by eliminating the redundant terms. To eliminate the redundant gene sets, we first surveyed all the MSigDB GOBP gene sets to identify which GOBP terms are contained in others, or which GOBP gene sets are subsets of other GOBP gene sets. After determining all the significant GOBP terms from the overrepresentation analysis, we removed the GOBP terms that are contained in other GOBP terms.
For the figures in the main text (Figures 2 and 3), we identified the top 3 GOBP terms for Protocol 2 for each SARS-CoV-2 proteins, because Protocol 2 produces the maximum network relative entropy due to SARS-CoV-2 PPIs. To estimate when a biological process may become a candidate for dysregulation due to SARS-CoV-2 PPIs, we considered the effect of finite responses and large deviations theory. We present this method for Protocol 1 results, where the reference state for each protein is {0.5, 0.5}, but this analysis can be performed for other reference states too [12,13]. As mentioned in SI: 4, we consider the proteins v i that have a relative entropy H ref (v i ) ≥ 0.01 bits to be part of the communicated set. The sum of the PPI and the communicated set is the total interaction set. A typical outcome of comparing the total interaction set against a particular BP gene set is shown in SI: 6.

Genes in BP
Genes not in BP Row total Total interaction set genes a = a PPI + a communicated b a + b Reference state genes  The p-value as a function of the variables in Table. SI: 6, which is the probability of obtaining an intersection greater than or equal to a by chance, is given by the hypergeometric survival function as p-value = S hg (a − 1, M, a + c, a + b) where S hg (·) is the hypergeometric survival function. If a protein is in the reference state {0.5, 0.5} and we sample or draw responses of the protein state for N times, then the empirical distribution {P (v ), may be greater than zero, without receiving any communication from SARS-CoV-2 PPIs. The probability of observing such a deviation can be determined using Sanov's theorem, which involves solving a constrained relative entropy minimization problem [13]. The solution for a reference state of {0.5, 0.5} is particularly simple: the probability of obtaining an empirical state after N responses that has a relative entropy greater than or equal to H ref (v i ) is 2 −H ref (vi)N . For example, NSP14 communicates to 5 proteins under Protocol 1 (Table. SI: 3). In the absence of SARS-CoV-2 PPIs, the probability of the relative entropy of the empirical state to be H ref (v i ) or more is in   To determine the significant biological processes as a function of N , we performed a Monte Carlo (MC) simulation to determine the N -dependent Q values for the GOBP terms identified in SI: 5 under Protocol 1. For each MC sample, we eliminated each protein v i from the communicated set with the probability 2 −H ref (vi)N , which means the total interaction set was itself a random variable corresponding to each MC sample. Subsequently, the outcome of comparing a sample total interaction set with a BP gene set (Table SI: 6) results in random variablesã,b,c andd (still satisfying the constraints mentioned in Table SI: 6). So, from n S MC samples we computed the average p-value asp -value = 1 n S n S k=1 S hg (ã (k) − 1, M, a + c,ã (k) +b (k) ) ( The average p-value from Eq. (6) was then used to identify the significant biological processes after Benjamini-Hochberg correction with a false-discovery rate of 1%. The Q values in Fig. 3 of the main text and in Fig. SI: 6 were computed using this procedure. Fig. 3(b) shows the Q-N curves for the proteins NSP12 and ORF9c, respectively. The plots for the other SARS-CoV-2 proteins are in Fig. SI: 6. We used n S = 10 4 for this calculation and verified for convergence using n s = 4 × 10 4 .