Simulation of the Dynamics of Primary Immunodeficiencies in B Cells

Primary immunodeficiencies (PIDs) are a group of over 300 hereditary, heterogeneous, and mainly rare disorders that affect the immune system. Various aspects of immune system and PID proteins and genes have been investigated and facilitate systems biological studies of effects of PIDs on B cell physiology and response. We reconstructed a B cell network model based on data for the core B cell receptor activation and response processes and performed semi-quantitative dynamic simulations for normal and B cell PID failure modes. The results for several knockout simulations correspond to previously reported molecular studies and reveal novel mechanisms for PIDs. The simulations for CD21, CD40, LYN, MS4A1, ORAI1, PLCG2, PTPRC, and STIM1 indicated profound changes to major transcription factor signaling and to the network. Significant effects were observed also in the BCL10, BLNK, BTK, loss-of-function CARD11, IKKB, MALT1, and NEMO, simulations whereas only minor effects were detected for PIDs that are caused by constitutively active proteins (PI3K, gain-of-function CARD11, KRAS, and NFKBIA). This study revealed the underlying dynamics of PID diseases, confirms previous observations, and identifies novel candidates for PID diagnostics and therapy.

Introduction

The human immunome covers the entirety of genes and proteins that are essential for innate and adaptive immunity (1). Many of these proteins are involved in extensive interaction networks. We have previously defined and characterized the essential immunome interactome, i.e., the totality of interactions in the immune system (2). These data can be used for many purposes including studies of the dynamics of the immune response in health and disease. The immunome interactome is not stable; it varies between cell types and even within them depending on the timing and localization of expressed and active proteins.

B cells produce antibodies. The membrane-bound antibody component of the B cell receptor (BCR) recognizes foreign antigens and triggers a cascade of signal transduction events that lead to the activation as well as nuclear transport of specific transcription factors (TFs) (3). In the nucleus, these TFs activate transcription of B cell proliferation and response genes. B cells can be divided into subpopulations based on the expression of surface markers. B cell maturation is a complex process that proceeds from hematopoietic stem cells via pro- and pre-cell stages to immature and mature B cells and finally to plasma or memory cells. The B cell subtypes play different and sometimes overlapping roles in immune response by virtue of their differing master regulator TFs (4) and surface receptors (5–7).

Primary immunodeficiencies (PIDs) are intrinsic diseases of the immune system, mostly rare and typically with heterogeneous phenotypes. Already over 300 PIDs have been described. Disease-causing variants have been collected to IDbases (8) and some other databases. Differential diagnosis of PIDs is often difficult because of overlapping signs and symptoms, sometimes only a genetic test can provide the definitive answer. There are some classification schemes for PIDs, including the one from the International Union of Immunological Societies (IUIS) expert committee for PIDs (9) and the network-based classification that clusters the diseases based on signs, symptoms and laboratory parameters (10).

The reconstruction of cellular networks in systems biology facilitates simulations where diseases are modeled as perturbations or alterations (11, 12). Since experimental studies are very tedious, they have been limited to small networks. Mathematical network simulations can offer insight into the dynamics of biomolecular interactions in cellular processes both in health and disease. Some protein–protein interaction (PPI) networks in B cells and their role in various diseases have been studied (13–15). These interactions can be disrupted by disease-causing variations. Quantitative dynamics studies are computationally intensive when the number of parameters is extensive. Further problems occur because of lack of kinetic parameters and reaction constants. Therefore, other approaches have been developed to study larger networks of up to hundreds of nodes by using qualitative and semi-quantitative dynamics (16–18).

Here, we employed the semi-quantitative method of normalized HillCube Boolean approach (19) to simulate the dynamics during the activation of B cells. Variants that cause PIDs were used to study the effects of PID perturbations. We conducted synchronous update simulations and validated them in silico. The simulations reproduced known trends due to variations in PIDs in essential signal transduction pathways during B cell development from pre- to mature B cell (20). Moreover, we found several novel proteins affected by PIDs and detected novel PID candidates.

The content of this article is the full paper of part of a published doctoral thesis (21).

Results

The Primary BCR Network

To investigate B cell signal transduction, we started by building the central networks based on data from literature. Then we constructed reaction equations for the primary BCR activation and response and converted them into a Boolean network model (Table S1 in Supplementary Material). These interactions were manually defined as Boolean equations using the sum-of-product form similar to our previous study on T-cell networks (22). Proteins were treated as Boolean variables in the hypergraph. Activating (+) and inhibiting (−) interactions were represented as edges with pointed or blunt arrowheads, respectively. Each edge originates from a source and points to a target and indicates the direction of signal transduction (Figure 1). Multiple edges pointing at the same node represent summed interactions and were represented by the OR operator, while multiple incoming edges that together activate or inhibit a protein were connected to the target using the AND gate. Totally 17 input nodes in the network did not have incoming links (Figure 1).

FIGURE 1

Figure 1. B cell activation Boolean network model. The network consists of 222 links and 144 nodes (including the AND Boolean operator), 17 of which are input nodes that have no links pointing to them. The Boolean network represents the B cell activation events. The rectangles for non-primary immunodeficiencies (PIDs) are white and for PID proteins are gray. Ellipses denote the AND gates. Input nodes are represented by hexagonal boxes. Activating links have pointed arrows while inhibiting links have blunt heads. Signal 1 represents B cell receptor (BCR) complex. Since the network focuses on BCR, its coreceptor (CD19/CD21/CD81), and costimulatory receptor, CD40 signaling events, only major events, e.g., for survival signaling and response have been fully considered.

We analyzed the structure of the network and the signaling paths between the initial events of the BCR activation and the late events of major TFs required for the expression of response genes. The BCR complex, its coreceptor, CD19/21/81, and the CD40 costimulatory receptor are involved in the initial events, while the TFs ELK1, AP1, NFAT, and NF-κB control the late events of B cell activation (23).

The BCR is activated when it binds to an antigen (signal 1). Another signal (signal 2) through the costimulatory receptor, via CD40, cytokines and CD19/CD21/CD81 complex is needed to elicit survival, and response (24, 25). The multiple paths from receptors to TFs guarantee a fail-safe and robust B cell activation (26–28). The sensitivity of gene activation is likely modulated by different routes (29).

We identified paths from signals 1 and 2 to the major response TFs including ATF2, CREB1, ELK1, ETS1, FOXO1, JUN, MEF2C, MTOR, NFKB1, and NFAT. For this purpose, we converted the Boolean network into an interaction network (Figure 2) to capture the dependencies, interactions, and thus, the paths through which signals are transduced. The entire interaction network is a single strongly connected component of 97 nodes interconnected by 168 links (Table S1 in Supplementary Material).

FIGURE 2

Figure 2. Boolean model transformed into its underlying interaction graph. The network consists of nodes and links derived from the Boolean network model without the AND operator. The interaction graph contains 107 nodes and 188 links and represents the underlying interaction network of the model. The nodes are as described in Figure 1. The network shows the paths through which signals from the receptors are channeled through the network to the transcription factors, which turn on the response genes.

Signaling loops indicate the dynamic nature of a network. We analyzed proteins that are essential for transducing signals between the receptor components (BCR and CD40LG) to downstream TFs. Proteins whose Boolean update equations are along most of the loops were considered essential. In total, we found 6,542 such loops of which the longest spans 26 nodes and the shortest 2 nodes. The mean length of the loops is 15 nodes. All the loops included PID proteins (Table S2 in Supplementary Material). SYK, PLCG2, LYN, PI3, PI3K, IP3, and CA are the nodes that are most frequently identified along the loops. These results show that the known PID proteins are central for the B cell pathways and thus harmful variants in these proteins would be detrimental for the signal flow.

In Silico Validation of Reconstructed Network and Identification of the Wild-Type Attractor

The engagement of the BCR, its coreceptor, CD19/CD21/CD81 complex and the CD40 stimulatory receptor, triggers a series of signal transduction cascades in B cells (3), which are captured in our reconstructed network (Figure 1; Table S1 in Supplementary Material). The signaling cascades activate numerous TFs, including ATF2 (30), CREB1 (31), ELK1 (31), ETS1 (32), FOXO1 (33), JUN (34, 35), MEF2C (36), NFKB1 (37), and NFAT (38). The reconstructed network is considered cogent when the major TFs and the signaling components that lead to their activation are turned on.

To ensure that our model reproduces B cell activation, we performed a large number of simulations by iteratively modifying the parameters and initial states of the input nodes, and ensuring that the network represented the main signaling events. We used normalized HillCube dynamic simulations (19) with signals 1 and 2 turned on and validated the simulations in silico. In addition, we performed simulations by turning signal 1 on and signal 2 off and vice versa. When signal 1 was turned off, ATF2 and JUN were turned on, while all other TFs were turned off, indicating the importance of the signal flow via this route. This effect is corroborated by previous experiments showing that when the CD40 receptor is stimulated by interaction with CD40LG the MAPK signal transduction pathways are activated. ATF2 and JUN are both activated by the MAPK signaling pathways (39). Although all other TFs were turned on when signal 2 was off, NFKB1, a very important survival signal for B cells, was turned off. This result is supported by experiments showing that the NF-κB pathway is activated by signal 1 (40). These simulations were accompanied by the tuning of the parameters such that the simulation results should comply with previously published experimental. We used the default values of the parameters for all equations except for those in Table 1.

TABLE 1

Table 1. Tuned parameters of nodes in the Odefy-simulated B cell network model.

We ran the normalized HillCube simulations until the networks reached their attractor states. The signaling events triggered by BCR activation were run until the network settled in an attractor after about 80 updates. The resulting attractors were in accordance with published experimental results (3, 41). For example, upon BCR stimulation, a cascade of phosphorylation leads to the activation of Src kinases (LYN, BLK, or FYN), SYK, BLINK BTK, PLCG2, and other proteins that form the early activation complex, which in turn leads to the calcium signaling (42). The activation of these proteins is captured in the basin of attraction of Figure 3. The simulations indicate also the activation and regulation of the major downstream TFs (ATF2, CREB1, ELK1, ETS1, FOXO1, JUN, MEF2C, NFKB1, and NFAT) when the signals 1 and 2 are turned on. In conclusion, the simulated network transmits signals to the TFs when both signals 1 and 2 are on, just as expected. Perturbations of key factors impair this flow of information and indicate disease-related processes due to defective signaling.

FIGURE 3

Figure 3. Attractor basin of the B cell network model’s normalized HillCube simulation. The basin of attractors of the B cell network model was simulated using the normalized HillCube algorithm. The horizontal axis denotes updates in time steps in arbitrary units. The simulation was run until the network reached a point attractor (the activating attractor) after about 80 time steps. At time points 110–210 the modulators of the B cell receptor signaling events (INPP5D, FCGR2B, PTPN6, PTPN11, PTPN12, PTEN, and HDAC4) were transiently activated. The network reached another attractor (the modulating attractor) at around time point 230. The simulations then continued until the activating attractor was reached again around time point 290.

PID Failure Analysis

We studied the effects of disease-causing variations on the long-term dynamics of B cells by perturbing PID proteins in the network model. Twenty-three PIDs are known to affect the proteins in the network including BCL10, BLNK, BTK, gain-of-function (GOF), and loss-of-function CARD11, CD19, CD21, CD40, CD81, IKKB, KRAS, LYN, MALT1, MS4A1, NEMO, NFKB1, NFKBIA, ORAI1, PI3K, PLCG2, PTPRC, STIM1, and WIPF1 deficiencies. Although no case of LYN deficiency has been reported, LYN was studied to represent PIDs connected to the Src-family kinases (SFKs). Thus, BLK was represented by LYN in the network model. To the best of our knowledge, BLK is the only SFK involved in B cell PIDs (43). The known PID proteins were identified from the ImmunoDeficiency Resource (44), IDbases (8), the classification by the IUIS expert committee for PIDs (9), and a recent review (45). These proteins are expressed at different stages during the B cell development. Here, we focused on PIDs that occur during pre- and mature B cell developmental stages. The effects of complete knockouts or knockins of these proteins were investigated by turning them off or on during simulations, respectively. The resulting perturbed attractors were compared with that of the wild type.

None of the major TF signaling pathways were disrupted in the CARD11, KRAS, and PI3K overexpression PID attractors, and CD19, CD81, and WIPF1 knockout attractors (Figure 4). Several TF pathways were dysregulated in the attractors of PIDs involved in the early events of the BCR-dependent B cell activation including BLNK, BTK, CD21, CD40, LYN, MS4A1, and PTPRC (Figure 4). The perturbations of most of the PIDs indicate profound effects in the pathways essential for B cell survival and response. In the wild-type attractor, the ETS1 TF is turned on during the simulation. This is in accordance with its inhibitory role in BCR response (46). Furthermore, ETS1-deficient mice show enhanced expression of activating markers and increased secretion of autoantibodies (47). However, ETS1 is turned off in the CD21, CD40, MS4A1, ORAI1, PLCG2, PTPRC, and STIM1 PID attractors. CREB1 is associated with antigen-BCR-dependent survival signals. The activating TFs CREB1, MEF2C, and NFAT pathways were dysregulated in the CD21, CD40, LYN, MS4A1, ORAI1, PLCG2, PTPRC, and STIM1 deficiency attractors. Except for CD19, CD18, and WIPF1 PID attractors, NF-κB signaling pathway was disrupted in all PID attractors. FOXO1, a critical downstream effector of the PI3K/AKT signal transduction axis involved in controlling growth arrest and apoptosis, was turned on in the wild-type simulation, as verified by experiments (48). In the PID simulations, the characteristics of FOXO1 are similar to that of the wild type except for CD21, CD40, LYN, MS4A1, and PTPRC PID attractors. ELK1, ATF2, and JUN pathways are partially dysregulated in the CD40 deficiency attractors, and normal in all other PID attractors. These TFs are phosphorylated and activated by the MAPK signaling relay and are associated with BCR-dependent survival signals (3).

FIGURE 4

Figure 4. Wild-type and primary immunodeficiency (PID) attractors of the B cell network simulation. The node states for the wild-type and the PID-perturbed attractors (knockout perturbation of BCL10, CARD11 loss-of-function, CD19, CD21, CD40, CD81, IKKB, LYN, MALT1, MS4A1, NEMO, NFKB1, NFKBIA, ORAI1, PI3K, PLCG2, STIM1, and WIPF1 and knockin perturbation of CARD11 gain-of-function, KRAS, NFKBIA, and PI3K) attractors. The attractors are represented by the rows while the states of the nodes in the attractors are represented on the columns. The state of a node for an attractor is represented by the color of the cell on the row of the attractor; black means inactive whereas white means active. The lighter the color of the cell, the more active the protein.

The simulation results agree with experimental evidence. We have chosen examples from the different parts of the network. CD40, variants which can lead to hyper-IgM syndrome (HIGM3) (49, 50) due to impaired somatic hypermutation generating antibody variation, provides one of the central inputs for B cell activation. In our simulations, CD40 defect has a profound effect on the entire network in line with its central role for BCR signaling. BTK is in the middle of the network. Variants in this protein lead to X-linked agammaglobulinemia due to defective BTK tyrosine kinase activity. Btk is critical for B-cell development, differentiation, and signaling. Harmful BTK variants block B-cell differentiation to pro- and pre-B-cell stages. BTK is involved in calcium mobilization and fluxes, cytoskeletal rearrangements, and transcriptional regulation involving NF-κB (51, 52). The simulation results indicate effects on NF-κB and IKK similar to experimental evidence.

NEMO is an inhibitor of κB kinase gamma and regulates NF-κB by phosphorylating IκB leading to degradation of it and subsequent activation of NF-κB (53). Thus, the results indicating effects on NFKB1 and IKK agree with experimental studies. Defective NEMO signaling causes X-linked hypohidrotic ectodermal dysplasia with immunodeficiency and other diseases including osteopetrosis and lymphedema.

Correlation With PID Severity

Primary immunodeficiencies differ greatly in severity from very mild to life-threatening conditions. Severe combined immunodeficiency (SCID) is associated with high vulnerability to infectious diseases and can be lethal (54). According to the IUIS classification (9), some of the PIDs in this study are associated with SCIDs with reduced numbers or absent T and B cells subtypes or isotypes. These include BCL10, CARD11, CD40, IKKB, MALT1, and PTPRC deficiencies. The attractors of these PIDs show less dysregulation except for CD40 (Figure 4). However, most of the PIDs disrupt the NF-κB pathways.

Combined immunodeficiencies (CIDs) are less severe compared with SCIDs and have variable clinical phenotypes. NEMO, GOF NFKBIA, ORAI1, STIM1, and WIPF1 constitute CIDs. Most of these proteins are components of calcium signaling, and impair, as expected, the NF-κB and NFAT pathways, both of which are activated by calcium signaling.

Gain-of-function variants in the PIK3CB gene that codes for the catalytic subunit of the PI3K heterodimeric complex are associated with a mild PID (55–57). In addition, variants in the gene that code for NFKBIA are associated with various forms of ectodermal dysplasia with immunodeficiency (58–61). The attractors for both PI3K and NFKBIA are similar to the wild-type indicative of a mild effect.

Many B cell PIDs are antibody deficiencies and have less severe but recurrent infections. These include BTK, BLNK, CARD11, CD19, CD21, CD81, PI3K, and STIM1 deficiencies. These PIDs show minor alterations to the investigated network, except for CD21 and STIM1. Nonetheless, the NF-κB pathway is disrupted in most of these PID simulations, as is the case with the experimental studies. Altogether, the PID severity analysis verifies that proteins involved in severe diseases are central for the network.

Novel PID Candidate Proteins

New PIDs are still being discovered, but due to their large number, rarity and overlapping symptoms, their diagnosis may be late, challenging and costly. Several classifications of PIDs have been introduced (9, 10), and candidate genes and proteins have been suggested (22, 62–65). Proteins that affect several pathways and are captured by the signaling loops are likely involved in PIDs. Proteins that appear in at least 10% of the loops include the majority of the investigated PIDs (15 of 22) and several proteins essential for B cell activation and functions. Interestingly, some of these proteins also indicate abrogated signaling in the attractors for most of the PIDs. To evaluate in silico the effects of perturbing proteins that have not, thus far, been experimentally connected to PIDs in Table 2, we performed knockout simulations for each of them. All the perturbed nodes impaired BCR-dependent survival and response signaling. Moreover, we analyzed from the Human Gene Connectome (HGC) (66) web server biologically and functionally close genes for known PIDs and found connections between proteins present in numerous loops and significant associations to known PID proteins (Table 2). Among these are SYK, PRKCB, IP3R, and GAB1 that are important for BCR activation and survival response. In addition, many of our candidate genes have been suggested in other studies (62, 63).

TABLE 2

Table 2. Novel primary immunodeficiency (PID) candidates.

Discussion

In this study, we reconstructed the B cell network model based on the literature, refined and in silico validated it, and then used it to study the semi-quantitative dynamic effects of PID knockouts and knockins. The model captures the main BCR activation signaling components in B cells. The normalized HillCube approach in the Odefy toolbox was used to simulate both the normal and the PID perturbations in the model for the B cell network. The perturbations qualitatively replicated experimental data for several PIDs at the pre- and mature B cell developmental stages.

The validity of the network and our approach were tested at several levels. First, the signal transduction from the BCR to TFs was shown to be intact and affected by perturbations of the proteins. Second, in many PIDs proteins were shown to be highly connected and involved in numerous signaling loops. Third, the resulting attractors from PID simulations were in line with experimental results for diseases in these proteins. Fourth, severity of several PIDs correlates with effects to information flow in the network. In summary, the network and simulations capture important characteristics of PIDs and can thus be used to extrapolate to other proteins in the network.

We compared the wild-type to the PID attractors and found severe signal transduction defects when CD21, CD40, LYN, MS4A1, ORAI1, PLCG2, PTPRC, and STIM1 were perturbed. The trends that indicate major changes in signaling patterns were captured in the knockout simulations. Minor differences were observed between the wild-type and CARD11, KRAS, and PI3K overexpression attractors, and CD19, CD81, and WIPF1 knockout attractors. The pathway for NF-κB was disrupted in all the PID knockouts except in the CD19, CD81, and WIPF1. These proteins connect receptor-dependent signals to the distal NF-κB pathway (67). knockout of any of these proteins may impair the IKK complex, the major NF-κB regulator, by leaving NFKBIA bound to NFKB1, thereby preventing its nuclear transportation and function as a transactivator (67).

After BCR-antigen activation, B cells undergo somatic hypermutation and antibody class switching. ELK1, ATF2, and JUN are effectors downstream to the RAS/MAPK pathways in the BCR signaling network. They control survival, differentiation and proliferation responses in B cells after activation (31). These TFs were not affected except (slightly) in the CD40 PID attractor. The ETS1 TF controls survival and differentiation, and thus its activity is increased through calcium ion signaling during BCR signaling (32, 46, 47, 68). In our simulations, the PID attractors for proteins involved in calcium signaling disrupted the ETS1 pathway. All other TFs that regulate survival and proliferation signals through the BCR signaling were abrogated in at least five PID attractors (Figure 4). These results verify that our simulation approach is effective when the affected proteins are at the core of the interconnected network or along non-redundant paths. No major changes were revealed in overexpression perturbations, as in PI3K, GOF CARD11, KRAS, and NFKBIA disorders. All the PID proteins, excluding those at input nodes, emerge in the signaling loops. LYN, STIM1, ORAI1, and PLCG2 perturbations, which cause major effects, are present in over two-thirds of the loops.

Our results show PID trends in the cellular dynamics of the B cells when the affected proteins are involved in non-redundant paths along the major TF signaling pathways. Perturbations of early events show more profound network disruption than those affecting late events. This is demonstrated in the profound outcomes when perturbing PID proteins taking part in the early events (CD21, CD40, CD81, LYN, MS4A1, and PTPRC), and less profound, but noticeable effects, in the intermediate and late events (IKKB, ORAI1, PLCG2, and STIM1).

The proteins in Table 2 that are not linked to known PIDs are essential for B cell activation and function and are affected in several B cell simulated attractors. Many of them have been previously identified as PID candidates. GAB1, GRB2, LAT2, MAP2K1, MAP3K7, MAPK14, MAPK3, MAPK8, PIK3AP1, PRKCB, RAF1, and SYK have been designated as candidate PID genes (63). Fourteen out of these 17 proteins were predicted as PID candidates in another study (62) including IKKA, GAB1, GRB2, LAT2, MAP2K1, MAP3K7, MAPK14, MAPK3, MAPK8, PIK3AP1, PRKCB, RAF1, and SYK. Several of the candidate genes are connected to PID proteins in the HGC (Table 2), providing independent proof for their significance.

According to the Panther databases (69), eight of the proteins in Table 2 are non-receptor protein kinases (GRB2, IKKA, MAP3K7, MAPK3, MAPK8, MAPK14, RAF1, and SYK), three are mitogen-activated and phosphotyrosine binding kinases (MAPK3, MAPK14, and RAF1), four have calcium ion binding activity (ITPR1, LAT2, MAP3K7, and PRKCB), and one has SH3/SH2 adaptor activity (GAB1). According to the Genecards and Malacards suites (70, 71), eight of the proteins are associated with various cancer types (GRB2, MAPK3, MAPK8, MAPK14, PRKCB, RASA1, SHC1, and SYK), four to different forms of cardiovascular diseases (MAP2K1, RAF1, RASA1, and SHC1), three to various Noonan syndromes (GAB1, MAP2K1, and RAF1), two to liver diseases (MAPK8 and SHC1), two to neuroendocrine neoplasm (MAPK3 and MAPK8), and one each to inflammatory bowel disease (SHC1), ectodermal dysplasia (IKKA), and Wiskott-Aldrich syndrome (GRB2). These proteins or their encoding genes could also have effects on PIDs. The listed proteins are strong PID candidates; however, their involvement in PIDs is yet to be experimentally verified.

Primary immunodeficiency candidate genes have been proposed in several studies (62–64). With reconstructed PPI network of immune system-specific proteins, proteins with high network statistics and PID-associated Gene Ontology term enrichment scores were considered as candidates (64). Itan and Casanova identified the top 1% of genes biologically close to known PIDs and selected those with similar Gene Ontology terms as the known PIDs (62). Support vector machine, a supervised machine learning technique, has been used to identify candidate PIDs (63). Several of the detected candidate genes above have been later confirmed to be PID-associated. Our approach focuses on B cell-specific network model, B cell intrinsic PIDs, the proteins that are connected to the PIDs and their effects on the BCR signaling in silico.

Differential diagnosis and treatment of PIDs is still challenging. Our approach provides novel insight to the mechanisms of PIDs in immune response signaling and presents new candidates for therapy and diagnosis. Experiments to validate the proposed PID candidates are still to be conducted. Furthermore, detailed quantitative models of the knockin perturbations would shed more light on their effects and require more detailed experimental data.

Materials and Methods

Network Reconstruction and Analysis

The B cell network model was reconstructed through extensive literature mining. We focused on the major components required for BCR/CD40-dependent B cell activation signaling events. Furthermore, we included only signaling events for undifferentiated B cells. Boolean equations were constructed based on literature evidence. The network is available in SBML format (Table S3 in Supplementary Material) and at the NDEx network provenance repository (72) with the UUID: 2554db2d-7533-11e8-a4bf-0ac135e8bacf. The reconstructed Boolean network model was used for wild-type and PID-perturbed simulations using Odefy, a Matlab toolbox (19), with the normalized HillCube method (73). NetDS, a Cytoscape plugin (74), was used for identifying signaling loops in the underlying interaction graph of the model. Data analysis was done with the R software, version 3.2.3 (75), and igraph, a library for network and graph analyses in R (76). Cytoscape, version 3.5.1 (77) was used for network visualization.

The analyses and simulation protocol were essentially similar to what we used to investigate T-cell PIDs (22). Briefly, in the Boolean model the nodes represent N signaling molecule or protein variables, X₁, X₁, …, X_N. The value of a variable is either 0 or 1 (78). Proteins, x_i, are influenced by their neighbors, R_i = {X₁, X₁, …, X_N}. The value of a protein at time t is updated at time t + 1 from the values of its neighbors, R_i, with update function B: {0, 1^N}. The discrete and synchronous update (78, 79) was performed according to the following equation:

$x_i\left(t+1\right)=B_i\left(x_{i1}\left(t\right),x_{i2}\left(t\right),\dots ,x_{i{N}_i}\left(t\right)\right)\text{\hspace{0.17em}ε\hspace{0.17em}}\{0,1\},i=1,2,\dots ,N.$

The ordinary differential equation (ODE) equivalent of the Boolean update functions, where x_i takes values [0, 1] was computed using

${\dot{x}}_i=\frac{1}{{\text{τ}}_i}\left({\overline{B}}_i({\overline{x}}_{i1},{\overline{x}}_{i2},\dots ,{\overline{x}}_{i{N}_i}\right)-{\overline{x}}_i),$

where ${\overline{B}}_i$ is a continuous homolog of the discrete function B_i, parameter τ_i is the life-time of the protein, and ${\overline{x}}_i$ describes its decay.

Odefy (19) transforms the discrete Boolean to the continuous system of ODEs and computes the solution of the system using the BooleCubes (73) as follows:

${\overline{B}}^I\left({\overline{x}}_1,{\overline{x}}_2,\dots ,{\overline{x}}_N\right)={\displaystyle \sum _{x_1=0}^1}{\displaystyle \sum _{x_1=0}^1}\dots {\displaystyle \sum _{x_N=0}^1}\left[B\left(x_1,x_2,\dots ,x_N\right).{\displaystyle \prod _{i=1}^N}\left(x_i{\overline{x}}_i+\left(1-x_i\right)\left(1-{\overline{x}}_i\right)\right)\right].$

${\overline{B}}^I$, the BooleCube, is obtained from interpolating B_i, the Boolean discrete update function. Since molecular interactions are not switch-like in behavior, the Hill function,

$\raisebox{1ex}{$f\left(\overline{x}\right)-{\overline{x}}^n$}\!\left/ \!\raisebox{-1ex}{$\left({\overline{x}}^n-{\overline{k}}^n\right)$}\right.,$

was applied to the BooleCube to obtain a sinusoidal function, the HillCube (73) as follows, ${\overline{B}}^H({\overline{x}}_1,\dots ,{\overline{x}}_N)={\overline{B}}^I(f_1({\overline{x}}_1),\dots ,f_N({\overline{x}}_N))$.

The parameter n (the Hill coefficient) represents the cooperativity between the protein interactions, while parameter k represents the value at which the activation is half-maximal.

To obtain perfect homologs of the Boolean update functions B_i, the HillCube functions were normalized to the unit interval to give the normalized HillCube (73) as follows, ${\overline{B}}^{Hn}({\overline{x}}_1,\dots ,{\overline{x}}_N)={\overline{B}}^I\left(\frac{f_1({\overline{x}}_1)}{f_1(1)},\dots ,\frac{f_N({\overline{x}}_N)}{f_N(1)}\right)$.

Basin of Attraction and Attractor Identification

The Odefy toolbox was used to simulate the qualitative dynamics of the network model (19). We used the normalized HillCube functions with synchronous updates. The default parameters for the normalized HillCube were n = 3, k = 0.5, and τ = 1. Except for the nodes in Table 1, the default parameter values were used. The variable n represents the Hill exponent of the Hill function and is used for converting the discrete Boolean update functions that take value 0 or 1 into their continuous BooleCube equivalents that have values [0, 1]. It captures the influence that the nodes of the same Boolean equation have on each other. k is a variable to control the continuous relaxation of the Boolean step function. It represents the value at half-maximal activation of a protein. τ is a decay parameter; for each signaling molecule, the higher its value, the slower the decay of the molecule.

The parameters in Table 1 were adjusted so that the wild-type attractor represents experimental data from the literature. The simulations were run until the network dynamics settled in the attractor that captures BCR activation.

Perturbation

Analysis of PID effects was performed for each protein encoded by a PID gene using the normalized HillCube simulations. For each perturbation, the node was converted to an input before assigning a state, either off or on, depending on the PID. This state was maintained until the simulation transitioned into the attractor. The parameter values used in the wild-type simulations were maintained for all the PID-perturbed simulations. The result of the simulation is the perturbed PID attractor.

Primary Immunodeficiencies

Primary immunodeficiency proteins expressed in B cells were retrieved from the IDbases (8, 44), the updated IUIS expert committee classification of PIDs (9), and a recent survey (45), and used for the PID failure mode simulations. The studied PIDs included BCL10, BLNK, BTK, CARD11, CD19, CD21, CD40, CD81, IKKB, KRAS, LYN, MALT1, MS4A1, NEMO, NFKB1, NFKBIA, ORAI1, PI3K, PLCG2, STIM1, and WIPF1 deficiencies.

Author Contributions

GT designed and developed the systems biology approach, interpreted data, performed data processing, interpreted the results, and wrote the manuscript. MV designed the research question, interpreted data, corrected the manuscript, and supervised the work. GT and MV have read and approved the final version of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This study is part of a published doctoral thesis (21). This study was funded by Vetenskapsrådet and the Alfred Österlunds stiftelse.

Supplementary Material

The Supplementary Material for this article can be found online at https://www.frontiersin.org/articles/10.3389/fimmu.2018.01785/full#supplementary-material.

Table S1. B cell activation Boolean network model update equations. The table lists Boolean equations of protein activation used in the network model and simulations.

Table S2. Signaling loops from the B cell interaction network of the Boolean network. The table consists of the number signaling loops that each node was a part of. All the loops included primary immunodeficiency proteins. All input nodes were excluded from the table. Signal1 (antigen) and CD40LG (Signal 2) were the source of all loops, and the rest of the input signals did not occur in any loop.

Table S3. SBML qual. The B cell network model in SBML qual format.

References

1. Ortutay C, Siermala M, Vihinen M. Molecular characterization of the immune system: emergence of proteins, processes, and domains. Immunogenetics (2007) 59(5):333–48. doi:10.1007/s00251-007-0191-0

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ortutay C, Vihinen M. Efficiency of the immunome protein interaction network increases during evolution. Immunome Res (2008) 4:4. doi:10.1186/1745-7580-4-4

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Dal Porto JM, Gauld SB, Merrell KT, Mills D, Pugh-Bernard AE, Cambier J. B cell antigen receptor signaling 101. Mol Immunol (2004) 41(6–7):599–613. doi:10.1016/j.molimm.2004.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Schmidlin H, Diehl SA, Blom B. New insights into the regulation of human B-cell differentiation. Trends Immunol (2009) 30(6):277–85. doi:10.1016/j.it.2009.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Klinman NR, Press JL. The B cell specificity repertoire: its relationship to definable subpopulations. Transplant Rev (1975) 24:41–83.

Google Scholar

6. Leandro MJ. B-cell subpopulations in humans and their differential susceptibility to depletion with anti-CD20 monoclonal antibodies. Arthritis Res Ther (2013) 15(Suppl 1):S3. doi:10.1186/ar3908

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Meyer-Bahlburg A, Rawlings DJ. Differential impact of toll-like receptor signaling on distinct B cell subpopulations. Front Biosci (Landmark Ed) (2012) 17:1499–516. doi:10.2741/4000

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Piirilä H, Väliaho J, Vihinen M. Immunodeficiency mutation databases (IDbases). Hum Mutat (2006) 27(12):1200–8. doi:10.1002/humu.20405

CrossRef Full Text | Google Scholar

9. Picard C, Al-Herz W, Bousfiha A, Casanova JL, Chatila T, Conley ME, et al. Primary immunodeficiency diseases: an update on the classification from the International Union of Immunological Societies Expert Committee for Primary Immunodeficiency 2015. J Clin Immunol (2015) 35(8):696–726. doi:10.1007/s10875-015-0201-1

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Samarghitean C, Ortutay C, Vihinen M. Systematic classification of primary immunodeficiencies based on clinical, pathological, and laboratory parameters. J Immunol (2009) 183(11):7569–75. doi:10.4049/jimmunol.0901837

PubMed Abstract | CrossRef Full Text | Google Scholar

11. del Sol A, Balling R, Hood L, Galas D. Diseases as network perturbations. Curr Opin Biotechnol (2010) 21(4):566–71. doi:10.1016/j.copbio.2010.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Goh KI, Cusick ME, Valle D, Childs B, Vidal M, Barabasi AL. The human disease network. Proc Natl Acad Sci U S A (2007) 104(21):8685–90. doi:10.1073/pnas.0701361104

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Basso K, Margolin AA, Stolovitzky G, Klein U, Dalla-Favera R, Califano A. Reverse engineering of regulatory networks in human B cells. Nat Genet (2005) 37(4):382–90. doi:10.1038/ng1532

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mendez A, Mendoza L. A network model to describe the terminal differentiation of B cells. PLoS Comput Biol (2016) 12(1):e1004696. doi:10.1371/journal.pcbi.1004696

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Sciammas R, Li Y, Warmflash A, Song Y, Dinner AR, Singh H. An incoherent regulatory network architecture that orchestrates B cell diversification in response to antigen signaling. Mol Syst Biol (2011) 7:495. doi:10.1038/msb.2011.25

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Li FT, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci U S A (2004) 101(14):4781–6. doi:10.1073/pnas.0305937101

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Thakar J, Saadatpour-Moghaddam A, Harvill ET, Albert R. Constraint-based network model of pathogen-immune system interactions. J R Soc Interface (2009) 6(36):599–612. doi:10.1098/rsif.2008.0363

PubMed Abstract | CrossRef Full Text | Google Scholar

18. von Dassow G, Meir E, Munro EM, Odell GM. The segment polarity network is a robust development module. Nature (2000) 406(6792):188–92. doi:10.1038/35018085

CrossRef Full Text | Google Scholar

19. Krumsiek J, Poelsterl S, Wittmann DM, Theis FJ. Odefy – from discrete to continuous models. BMC Bioinformatics (2010) 11:233. doi:10.1186/1471-2105-11-233

CrossRef Full Text | Google Scholar

20. Hardy RR, Hayakawa K. B cell development pathways. Annu Rev Immunol (2001) 19:595–621. doi:10.1146/annurev.immunol.19.1.595

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Teku NG. Computational Analysis on the Effects of Variations in T and B Cells: Primary Immunodeficiencies and Cancer Neoepitopes. Ph.D. thesis. Lund, Sweden: Lund University (2017).

Google Scholar

22. Teku NG, Vihinen M. Simulation of the dynamics of primary immunodeficiencies in CD4+ T-cells. PLoS One (2017) 12(4):e0176500. doi:10.1371/journal.pone.0176500

PubMed Abstract | CrossRef Full Text | Google Scholar

23. LeBien TW, Tedder TF. B lymphocytes: how they develop and function. Blood (2008) 112(5):1570–80. doi:10.1182/blood-2008-02-078071

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Cherukuri A, Shoham T, Sohn HW, Levy S, Brooks S, Carter R, et al. The tetraspanin CD81 is necessary for partitioning of coligated CD19/CD21-B cell antigen receptor complexes into signaling-active lipid rafts. J Immunol (2004) 172(1):370–80. doi:10.4049/jimmunol.172.1.370

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Sato S, Miller AS, Howard MC, Tedder TF. Regulation of B lymphocyte development and activation by the CD19/CD21/CD81/Leu 13 complex requires the cytoplasmic domain of CD19. J Immunol (1997) 159(7):3278–87.

PubMed Abstract | Google Scholar

26. Rahman MM, McFadden G. Myxoma virus lacking the pyrin-like protein M013 Is sensed in human myeloid cells by both NLRP3 and multiple toll-like receptors, which independently activate the inflammasome and NF-κB innate response pathways. J Virol (2011) 85(23):12505–17. doi:10.1128/JVI.00410-11

CrossRef Full Text | Google Scholar

27. Tibrewal N, Wu Y, D’mello V, Akakura R, George TC, Varnum B, et al. Autophosphorylation docking site tyr-867 in mer receptor tyrosine kinase allows for dissociation of multiple signaling pathways for phagocytosis of apoptotic cells and down-modulation of lipopolysaccharide-inducible NF-κB transcriptional activation. J Biol Chem (2008) 283(6):3618–27. doi:10.1074/jbc.M706906200

CrossRef Full Text | Google Scholar

28. Xu WF, Nair JS, Malhotra A, Zhang JJ. B cell antigen receptor signaling enhances IFNγ-induced Statl target gene expression through calcium mobilization and activation of multiple serine kinase pathways. J Interferon Cytokine Res (2005) 25(2):113–24. doi:10.1089/jir.2005.25.113

CrossRef Full Text | Google Scholar

29. Guo X, Wang XF. Signaling cross-talk between TGF-β/BMP and other pathways. Cell Res (2009) 19(1):71–88. doi:10.1038/cr.2008.302

CrossRef Full Text | Google Scholar

30. Gupta S, Campbell D, Derijard B, Davis RJ. Transcription factor ATF2 regulation by the JNK signal transduction pathway. Science (1995) 267(5196):389–93. doi:10.1126/science.7824938

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yasuda T, Sanjo H, Pages G, Kawano Y, Karasuyama H, Pouyssegur J, et al. Erk kinases link pre-B cell receptor signaling to transcriptional events required for early B cell expansion. Immunity (2008) 28(4):499–508. doi:10.1016/j.immuni.2008.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Sementchenko VI, Watson DK. Ets target genes: past, present and future. Oncogene (2000) 19(55):6533–48. doi:10.1038/sj.onc.1204034

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hinman RM, Bushanam JN, Nichols WA, Satterthwaite AB. B cell receptor signaling down-regulates forkhead box transcription factor class O 1 mRNA expression via phosphatidylinositol 3-kinase and Bruton’s tyrosine kinase. J Immunol (2007) 178(2):740–7. doi:10.4049/jimmunol.178.2.740

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Davis RJ. MAPKs: new JNK expands the group. Trends Biochem Sci (1994) 19(11):470–3. doi:10.1016/0968-0004(94)90132-5

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Karin M, Hunter T. Transcriptional control by protein phosphorylation: signal transmission from the cell surface to the nucleus. Curr Biol (1995) 5(7):747–57. doi:10.1016/S0960-9822(95)00151-5

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wilker PR, Kohyama M, Sandau MM, Albring JC, Nakagawa O, Schwarz JJ, et al. Transcription factor Mef2c is required for B cell proliferation and survival after antigen receptor stimulation. Nat Immunol (2008) 9(6):603–12. doi:10.1038/ni.1609

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Francis DA, Karras JG, Ke XY, Sen R, Rothstein TL. Induction of the transcription factors NF-κb, Ap-1 and NFAT during B-cell stimulation through the CD40 receptor. Int Immunol (1995) 7(2):151–61. doi:10.1093/intimm/7.2.151

CrossRef Full Text | Google Scholar

38. Hao S, Kurosaki T, August A. Differential regulation of NFAT and SRF by the B cell receptor via a PLCγ-Calcium dependent pathway. EMBO J (2003) 22(16):4166–77. doi:10.1093/emboj/cdg401

CrossRef Full Text | Google Scholar

39. Sutherland CL, Heath AW, Pelech SL, Young PR, Gold MR. Differential activation of the ERK, JNK, and p38 mitogen-activated protein kinases by CD40 and the B cell antigen receptor. J Immunol (1996) 157(8):3381–90.

PubMed Abstract | Google Scholar

40. Dolmetsch RE, Lewis RS, Goodnow CC, Healy JI. Differential activation of transcription factors induced by Ca2+ response amplitude and duration. Nature (1997) 386(6627):855–8. doi:10.1038/386855a0

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Harnett MM, Katz E, Ford CA. Differential signalling during B-cell maturation. Immunol Lett (2005) 98(1):33–44. doi:10.1016/j.imlet.2004.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Baba Y, Kurosaki T. Role of calcium signaling in B cell activation and biology. In: Kurosaki T, Wienands J, editors. B Cell Receptor Signaling. Cham: Springer International Publishing (2016). p. 143–74.

Google Scholar

43. Compeer EB, Janssen W, van Royen-Kerkhof A, van Gijn M, van Montfrans JM, Boes M. Dysfunctional BLK in common variable immunodeficiency perturbs B-cell proliferation and ability to elicit antigen-specific CD4+ T-cell help. Oncotarget (2015) 6(13):10759–71. doi:10.18632/oncotarget.3577

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Samarghitean C, Väliaho J, Vihinen M. IDR knowledge base for primary immunodeficiencies. Immunome Res (2007) 3:6. doi:10.1186/1745-7580-3-6

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Vihinen M. Immunodeficiency, primary: affecting the adaptive immune system. eLS. Chichester: John Wiley & Sons Ltd (2015). p. 1–6.

Google Scholar

46. Luo W, Mayeux J, Gutierrez T, Russell L, Getahun A, Muller J, et al. A balance between B cell receptor and inhibitory receptor signaling controls plasma cell differentiation by maintaining optimal Ets1 levels. J Immunol (2014) 193(2):909–20. doi:10.4049/jimmunol.1400666

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Garrett-Sinha LA. Review of Ets1 structure, function, and roles in immunity. Cell Mol Life Sci (2013) 70(18):3375–90. doi:10.1007/s00018-012-1243-7

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Szydlowski M, Jablonska E, Juszczynski P. FOXO1 transcription factor: a critical effector of the PI3K-AKT axis in B-cell development. Int Rev Immunol (2014) 33(2):146–57. doi:10.3109/08830185.2014.885022

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Al-Saud BK, Al-Sum Z, Alassiri H, Al-Ghonaium A, Al-Muhsen S, Al-Dhekri H, et al. Clinical, immunological, and molecular characterization of hyper-IgM syndrome due to CD40 deficiency in eleven patients. J Clin Immunol (2013) 33(8):1325–35. doi:10.1007/s10875-013-9951-9

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Ferrari S, Giliani S, Insalaco A, Al-Ghonaium A, Soresina AR, Loubser M, et al. Mutations of CD40 gene cause an autosomal recessive form of immunodeficiency with hyper IgM. Proc Natl Acad Sci U S A (2001) 98(22):12614–9. doi:10.1073/pnas.221456898

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Mohamed AJ, Yu L, Backesjo CM, Vargas L, Faryal R, Aints A, et al. Bruton’s tyrosine kinase (Btk): function, regulation, and transformation with special emphasis on the PH domain. Immunol Rev (2009) 228(1):58–73. doi:10.1111/j.1600-065X.2008.00741.x

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Pal Singh S, Dammeijer F, Hendriks RW. Role of Bruton’s tyrosine kinase in B cells and malignancies. Mol Cancer (2018) 17(1):57. doi:10.1186/s12943-018-0779-z

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Shifera AS. The zinc finger domain of IKKgamma (NEMO) protein in health and disease. J Cell Mol Med (2010) 14(10):2404–14. doi:10.1111/j.1582-4934.2010.01054.x

CrossRef Full Text | Google Scholar

54. van der Burg M, Gennery AR. The expanding clinical and immunological spectrum of severe combined immunodeficiency. Eur J Pediatr (2011) 170(5):561–71. doi:10.1007/s00431-011-1452-3

CrossRef Full Text | Google Scholar

55. Angulo I, Vadas O, Garcon F, Banham-Hall E, Plagnol V, Leahy TR, et al. Phosphoinositide 3-kinase delta gene mutation predisposes to respiratory infection and airway damage. Science (2013) 342(6160):866–71. doi:10.1126/science.1243292

CrossRef Full Text | Google Scholar

56. Crank MC, Grossman JK, Moir S, Pittaluga S, Buckner CM, Kardava L, et al. Mutations in PIK3CD can cause hyper IgM syndrome (HIGM) associated with increased cancer susceptibility. J Clin Immunol (2014) 34(3):272–6. doi:10.1007/s10875-014-0012-9

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lucas CL, Kuehn HS, Zhao F, Niemela JE, Deenick EK, Palendira U, et al. Dominant-activating germline mutations in the gene encoding the PI(3)K catalytic subunit p110δ result in T cell senescence and human immunodeficiency. Nat Immunol (2014) 15(1):88–97. doi:10.1038/ni.2771

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Courtois G, Smahi A, Reichenbach J, Doffinger R, Cancrini C, Bonnet M, et al. A hypermorphic IκBα mutation is associated with autosomal dominant anhidrotic ectodermal dysplasia and T cell immunodeficiency. J Clin Invest (2003) 112(7):1108–15. doi:10.1172/JCI18714

CrossRef Full Text | Google Scholar

59. Janssen R, van Wengen A, Hoeve MA, ten Dam M, van der Burg M, van Dongen J, et al. The same IκBα mutation in two related individuals leads to completely different clinical syndromes. J Exp Med (2004) 200(5):559–68. doi:10.1084/jem.20040773

CrossRef Full Text | Google Scholar

60. Lopez-Granados E, Keenan JE, Kinney MC, Leo H, Jain N, Ma CA, et al. A novel mutation in NFKBIA/IKBA results in a degradation-resistant N-truncated protein and is associated with ectodermal dysplasia with immunodeficiency. Hum Mutat (2008) 29(6):861–8. doi:10.1002/humu.20740

PubMed Abstract | CrossRef Full Text | Google Scholar

61. McDonald DR, Mooster JL, Reddy M, Bawle E, Secord E, Geha RS. Heterozygous N-terminal deletion of IκBα results in functional nuclear factor κB haploinsufficiency, ectodermal dysplasia, and immune deficiency. J Allergy ClinImmunol (2007) 120(4):900–7. doi:10.1016/j.jaci.2007.08.035

CrossRef Full Text | Google Scholar

62. Itan Y, Casanova JL. Novel primary immunodeficiency candidate genes predicted by the human gene connectome. Front Immunol (2015) 6:142. doi:10.3389/fimmu.2015.00142

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Keerthikumar S, Bhadra S, Kandasamy K, Raju R, Ramachandra YL, Bhattacharyya C, et al. Prediction of candidate primary immunodeficiency disease genes using a support vector machine learning approach. DNA Res (2009) 16(6):345–51. doi:10.1093/dnares/dsp019

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Ortutay C, Vihinen M. Identification of candidate disease genes by integrating gene ontologies and protein-interaction networks: case study of primary immunodeficiencies. Nucleic Acids Res (2009) 37(2):622–8. doi:10.1093/nar/gkn982

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Piro RM, Di Cunto F. Computational approaches to disease-gene prediction: rationale, classification and successes. FEBS J (2012) 279(5):678–96. doi:10.1111/j.1742-4658.2012.08471.x

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Itan Y, Zhang SY, Vogt G, Abhyankar A, Herman M, Nitschke P, et al. The human gene connectome as a map of short cuts for morbid allele discovery. Proc Natl Acad Sci U S A (2013) 110(14):5558–63. doi:10.1073/pnas.1218167110

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Mitchell S, Vargas J, Hoffmann A. Signaling via the NF-κB system. Wiley Interdiscip Rev Syst Biol Med (2016) 8(3):227–41. doi:10.1002/wsbm.1331

CrossRef Full Text | Google Scholar

68. Russell L, John S, Cullen J, Luo W, Shlomchik MJ, Garrett-Sinha LA. Requirement for transcription factor Ets1 in B cell tolerance to self-antigens. J Immunol (2015) 195(8):3574–83. doi:10.4049/jimmunol.1500776

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Mi H, Lazareva-Ulitsky B, Loo R, Kejariwal A, Vandergriff J, Rabkin S, et al. The PANTHER database of protein families, subfamilies, functions and pathways. Nucleic Acids Res (2005) 33(Database issue):D284–8. doi:10.1093/nar/gki078

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Rappaport N, Twik M, Plaschkes I, Nudel R, Iny Stein T, Levitt J, et al. MalaCards: an amalgamated human disease compendium with diverse clinical and genetic annotation and structured search. Nucleic Acids Res (2017) 45(D1):D877–87. doi:10.1093/nar/gkw1012

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, et al. The GeneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinformatics (2016) 54 1.30.1–1.30.33. doi:10.1002/cpbi.5

CrossRef Full Text | Google Scholar

72. Pratt D, Chen J, Welker D, Rivas R, Pillich R, Rynkov V, et al. NDEx, the network data exchange. Cell Syst (2015) 1(4):302–5. doi:10.1016/j.cels.2015.10.001

CrossRef Full Text | Google Scholar

73. Wittmann DM, Krumsiek J, Saez-Rodriguez J, Lauffenburger DA, Klamt S, Theis FJ. Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling. BMC Syst Biol (2009) 3:98. doi:10.1186/1752-0509-3-98

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Le DH, Kwon YK. NetDS: a Cytoscape plugin to analyze the robustness of dynamics and feedforward/feedback loop structures of biological networks. Bioinformatics (2011) 27(19):2767–8. doi:10.1093/bioinformatics/btr466

PubMed Abstract | CrossRef Full Text | Google Scholar

75. R-Core-Team. R: A Language and Environment for Statistical Computing. (2016). Available from: http://www.R-project.org

Google Scholar

76. Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal, Complex Systems (2006) 1695(5):1–9.

Google Scholar

77. Kohl M, Wiese S, Warscheid B. Cytoscape: software for visualization and analysis of biological networks. Methods Mol Biol (2011) 696:291–303. doi:10.1007/978-1-60761-987-1_18

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Boole G. The calculus of logic. Cambridge Dublin Math J (1848) 3:183–98.

Google Scholar

79. Thomas R. Boolean formalization of genetic control circuits. J Theor Biol (1973) 42(3):563–85. doi:10.1016/0022-5193(73)90247-6

CrossRef Full Text | Google Scholar