From Discrete to Continuous Modeling of Lymphocyte Development and Plasticity in Chronic Diseases

The molecular events leading to differentiation, development, and plasticity of lymphoid cells have been subject of intense research due to their key roles in multiple pathologies, such as lymphoproliferative disorders, tumor growth maintenance and chronic diseases. The emergent roles of lymphoid cells and the use of high-throughput technologies have led to an extensive accumulation of experimental data allowing the reconstruction of gene regulatory networks (GRN) by integrating biochemical signals provided by the microenvironment with transcriptional modules of lineage-specific genes. Computational modeling of GRN has been useful for the identification of molecular switches involved in lymphoid specification, prediction of microenvironment-dependent cell plasticity, and analyses of signaling events occurring downstream the activation of antigen recognition receptors. Among most common modeling strategies to analyze the dynamical behavior of GRN, discrete dynamic models are widely used for their capacity to capture molecular interactions when a limited knowledge of kinetic parameters is present. However, they are less powerful when modeling complex systems sensitive to biochemical gradients. To compensate it, discrete models may be transformed into regulatory networks that includes state variables and parameters varying within a continuous range. This approach is based on a system of differential equations dynamics with regulatory interactions described by fuzzy logic propositions. Here, we discuss the applicability of this method on modeling of development and plasticity processes of adaptive lymphocytes, and its potential implications in the study of pathological landscapes associated to chronic diseases.

1. Introduction

The extensive accumulation of data from short and large-scale experiments involving a wide spectrum of biological functions of B and T lymphocytes in both, normal and pathological scenarios, has inspired an intensive research on molecular events leading to their early development, plasticity and emergency differentiation. As a result, the construction of regulatory networks has become a resourceful tool for the systems-level analyses of cell fate decisions through interconnection of molecular elements, such as biochemical signals provided by the microenvironment (e.g., cytokines, growth factors, transmembrane ligands, antigens, etc.) and transcriptional modules underlying the regulation of lineage-specific gene expression. Getting insights into the dynamical behavior of regulatory networks in biology requires simulation as continuous or discrete models (1). Discrete modeling, represented by Boolean and multi-valued network models, has been useful in differentiation processes of adaptive B and T lymphocytes (2–8), for molecular switching in cellular specification (9), for the prediction of microenvironment-dependent cell plasticity (6, 10), and for the analyses of signaling events occurring downstream activation of antigen recognition receptors (11, 12). Moreover, Boolean algebra has been used in cytometry to create combined gates for the identification and selection of cellular subsets and lymphoid phenotyping (13). Nevertheless, the utility of discrete models is limited as they cannot predict outcomes from quantitative biological experiments when working on phenomena sensitive to graded expression of transcription factors or biochemical gradients. This is the case of most diseases where lymphocytes are involved and non-discrete fluctuations in the microenvironment may influence cell differentiation and plasticity, affecting immune responses at the progression of chronic pathologies, such as lymphoproliferative disorders, tumor growth, diabetes, cardiovascular, and chronic respiratory diseases, among others. Discrete models might be then transformed into differential equations to allow a dynamical analyses of regulatory networks, as transformed continuous models, with potential implications in lymphoid cell- associated pathologies (14–17).

Here we propose the fuzzy logic transformation of a discrete model into a continuous model to compensate their disadvantages and to simulate biological systems with a well-known network architecture strongly influenced by concentration-dependent cues (Table 1).

TABLE 1

Table 1. Mathematical dynamic modeling subtypes: advantages and disadvantages.

2. Discrete Modeling of Lymphoid Differentiation Landscape

2.1. Boolean Interpretation of Molecular Data

To deeply understand the gene regulatory processes involved in cellular development, C. H. Waddington introduced in 1957 the metaphoric concept of epigenetic landscape (18). He proposed a unique perspective of cellular development as a ball rolling down within a landscape formed by peaks and valleys. Following its trajectory, the ball may finally fall into a valley, representing its final position that defines a steady-state -and a cellular fate-, also known as attractor. Waddington's epigenetic landscape was formalized, among others, by S. A. Kauffman, who studied the behavior of large networks of randomly interconnected binary “genes” with a dichotomous (on-off) behavior, establishing the principles of Boolean modeling (19). The assumption of a discrete transcriptional regulation was further investigated in Drosophila embryogenesis, showing that the gradient of Bicoid morphogen resulted from averaging binary states of transcriptional activity, active or inactive, at individual nuclei level (20).

The general system's behavior and the number of attractors of a Boolean or multi-valued regulatory network depends on topological characteristics, such as the number of components and the degree of interconnectivity among them. It is now recognized that biological networks are scale-free systems, which means that the nodes have a high diversity of number of edges, including few elements with many links and many elements with few links (21, 22). Scale-freeness provides, among other attributes: network robustness, better information spreading performance, and the property that the number of attractors is almost independent from the number of nodes (23, 24).

Mathematical modeling based on Boolean regulatory networks (BRN) provides meaningful qualitative information on the basic topology of relations that determine alternative cell fates and may be used for the analysis of biological circuits without requiring explicit values of the network parameters. In this type of approach, the network nodes represent genes, transcription factors, proteins mediating signaling cascades, RNA, environmental factors, etc., and links representing positive or negative regulation between pairs of nodes. The state variable of each node takes a discrete value of 0 (inhibited, or inactive) or 1 (expressed or active) (1). The state of each node at time t + 1 is specified by a dynamic mapping that depends on the state of its regulators at a previous time t:

$\begin{array}{ll}{q}_k(t+1)=F_k\left(q_1(t),\dots ,q_n(t)\right)& (1)\end{array}$

where F_k is a discrete function representing a logical proposition, also known as Boolean rules, constituted by elementary terms related by the logical connectives: AND (∧), OR (∨), and NOT (¬). Logical propositions satisfy Boolean's axiomatics, which complies associativity, commutativity, distributivity, absorptivity, and identity. The discrete nature of the truth values involved in Boolean logic propositions implies that this approximation is not always enough to investigate the enormous variability inherent to biological processes.

The dynamics induced by the Boolean mapping is completely determined once a set of initial expression values of the network components is specified. From a given initial set, the network nodes iteratively update their value based on the Boolean transfer rules until eventually reaching a steady-state determined by condition q_k(t + 1) = q_k(t). This latter condition specifies a fixed-point attractor. Then, the dynamics of a model is evaluated by tracking the trajectories from all the possible initial configurations in the states space toward the attractors. The size of the states space of a model is given by Ω= 2ⁿ where n is the number of nodes in the network. Alternatively, a cyclic attractor associated to the condition q_k(t + N) = q_k(t) may also arise after the simulation of some regulatory networks, where the integer number N signals the period of the attractor. Cyclic attractors are generally interpreted as oscillatory behaviors and are sustained by at least one negative feedback circuit in the network topology, which involves an odd number of inhibitory interactions (25). This type of attractors can be directly associated to biological events, for example, in models predicting cell cycle oscillations (26–28) or, sometimes they can be interpreted as intermediate or oscillatory activations in multi-valued and Boolean differentiation models, as has been reported with T cell attractors (7, 29). Each fixed-point and cyclic attractor is reached from a number ω of different initial conditions. The parameter ω denotes the size of the attraction basin which may be visualized as a ratio of areas in the epigenetic landscape. Consequently, the probability that a steady state is expressed is given by p = ω/Ω.

To briefly exemplify how a Boolean model is constructed we used the information compilated by Bhattacharya et al. (30) of the transcriptional core orchestrating the terminal differentiation of B cells into antibody-secreting plasma cells upon antigenic stimulation. The transcription factors to be considered were Pax-5, Bcl-6, and Blimp-1. Construction of the gene regulatory network and the Boolean transfer rules are based on evidence showing the existence of a mutual repression by Bcl-6 and Blimp-1, as do Blimp-1 and Pax-5, establishing a system with two double-negative feedback loops. Pax-5 and Bcl-6 are two transcriptional factors of high expression in B cells, down-regulated by Blimp-1 after its AP-1 mediated activation. In turn, AP-1 is phosphorylated downstream B cell stimulation with lipopolysaccharides. Beside the direct inhibition of Blimp-1, Bcl-6 can also act as a passive repressor through its binding to AP-1, blocking its transcriptional activity (Figure 1A). Such information is sufficient to predict two fixed-point attractors interpretable as B-cell and plasma cell configurations. The presence of at least one positive loop containing an even number of inhibitory regulations is necessary for the generation of multiple steady states (25). This type of models has been useful to merge independently published data from different molecular circuits involved in cellular specification, to probe how these circuits orchestrates differentiation, and to generate new testable hypothesis on missing interactions or cellular transitions.

FIGURE 1

Figure 1. (A) Boolean modeling of the transcriptional core regulating naive B cell to plasma cell (PC) differentiation. Inhibitory and activation interactions are represented in the network with truncated red lines and green arrows, respectively. (B) Epigenetic landscape of hematopoietic differentiation where valleys represent stable cellular states, however other cellular phenotypes may be represented as transitory stages. HSC, hematopoietic stem cell; MPP, multipotent progenitor; CMP, common myeloid progenitor; LMPP, lymphoid-primed multipotent progenitor; MEP, megakaryocyte/erythroid progenitor; MegP, unipotent megakaryocyte progenitor; GMP, granulocyte/macrophage progenitor; CLP, common lymphoid progenitor; NK/ILC, natural killer/innate lymphoid cell; ALP, all-lymphoid progenitors; BLP, biased lymphoid progenitors; ETP, early thymic progenitor; DN, double negative; DP, double-positive.

2.2. Genetic Regulatory Networks Underlying Lymphoid Specification

As of the discovery of HSCs by Ernest A. McCulloch and James E. Till in the 1960s, the hematopoietic system has served as the most recurrent biological model for the study of stem cell biology and differentiation. For many years, the differentiation process was represented as a hierarchical dichotomic model of strict myeloid/lymphoid branching. However, multiple observations mostly based on single cell experiments have challenged this classical view and introduced cell differentiation as a process of continuous transitions directed by two events running in parallel: the gradual commitment through the acquisition of lineage-specific features and the gradual lost of potential to generate cells of a different lineage (31–36) (Figure 1B).

In the metaphorical Waddington's view, the cell type positioned in the “summit” of the hematopoietic epigenetic landscape is the hematopoietic stem cell (HSCs) population, which resides in specialized niches within the bone marrow. Early specification begins upon “ball rolling” from the HSCs to the multipotent progenitor (MPP) attractors, either committing to myeloid or lymphoid lineages by differentiating into common myeloid progenitors (CMPs) or lymphoid-primed multipotent progenitors (LMPPs), respectively (37–39). As more is deciphered on the transcriptional network underlying the lymphoid differentiation, more is discovered about intermediate steps and novel transitional cell subpopulations. It is now well-known that LMPPs contain a mixture of myeloid and lymphoid-restricted progenitors, including early lymphoid progenitors (ELPs), giving rise to common lymphoid progenitors (CLPs), endowed with the ability of generating all types of adaptive and innate lymphocytes without noticeable myeloid potential, and some categories of dendritic cells (DCs) (40–50). The CLP population bisects into all-lymphoid progenitors (ALPs) and B-cell-biased lymphoid progenitors (BLPs) (44) that predominantly generate T and B lymphocyte precursors, respectively. From the ALP pool, some circulating progenitors reach the thymus and differentiate into early thymic progenitors (ETPs), progress to CD4/CD8 double-negative 2 cells (DN2) and DN3 cells. CD4/CD8 double-positive (DP) cells are then produced before differentiation toward CD4 or CD8 single-positive (SP) T effector cells (51). B cells reach also a partial maturation in the bone marrow (BM), following a series of sequential differentiation steps from prepro-B, pro-B, early pre-B and pre-B stages, where the rearrangement of immunoglobulin heavy-chain (IgH) genes takes place and results in the expression of the pre-B-cell receptor (pre-BCR). Downstream the pre-BCR activation and the signaling cascade deriving in a clonal expansion and the subsequent cell cycle arrest, a second wave of recombinases Rag1 and Rag2 expression induces the rearrangement of the immunoglobulin light-chain (IgL), marking the transition from pre-B-cell to immature B cells (44, 52, 53). Upon migration to the secondary lymphoid organs, T and B lymphocytes are exposed to antigens and signals provided by a number of immune cells in the microenvironment.

Even though differentiation transitions are now recognized as continuous processes, commitment to stable phenotypes is dependent on molecular switches that act as lineage-determining steps, what has made the differentiation process a target for its simulation through discrete models. More specifically, hematopoietic differentiation has been approached with discrete models at many levels (Figure 2), from the top of the epigenetic landscape hill sloping down to the final stages of mature cells production. The main type of information provided by the construction and simulation of BRN is obtained after the confirmation of the functional integration of the proposed components. This generally occurs validating the attractors and transitions with previous experimental observations. To compare with experimental data, Boolean models are subjected to different types of perturbations including permanent mutations (e.g., gene knock-out or overexpression), or temporal changes in the nodes activation value which can be understood as triggering cues for network state transitions. An example of this type of evaluations is the case of the hematopoietic stem/progenitor (HSPC) network model generated by Bonzanni. The HSPC model contains ten genes expressed in the immature stem cell population besides GATA1, which is expressed in the early progenitor MPP (2). The dynamic simulation of the HSPC network generated two single state attractors, one with an erythroid cell profile, and one with a non-hematopoietic cell profile with all genes turned-off, as well as a periodic attractor composed of 32 interconnected states with oscillatory activation values for four genes (Gata2, Zfpm1, Erg, and Eto2a) compatible with single cell gene expression data from HSPCs (54). The activation state of one or more genes in the states comprising the HSPC complex attractor were modified to compute the dynamic transitions and mapping the developmental route from HSPC toward erythrocyte, granulocyte, monocyte, natural killer (NK), B cell, CD4, or CD8 T cells profiles (55). This type of evaluation provided information about the stability of the HSPC attractor, the type of genes involved in the developmental route considering those that trigger differentiation, and the suggestion that there were missing interactions or components that avoid differentiation reversal.

FIGURE 2

Figure 2. Cellular attractors and transitions from the hematopoietic landscape reproduced through discrete modeling. Each color represents an independent article of hematopoietic discrete modeling, black arrows represent the direction of hematopoiesis toward the myeloid and lymphoid lineage. White nodes represent predicted phenotypes that have not been associated to experimental findings.

Furthermore, the myeloid/lymphoid branching has been addressed through the assembly of a GRN integrating 23 nodes that, when computed using a logical multi-valued formalism, produced four stable stages corresponding to CLPs, B-lineage cells, granulocyte-monocyte progenitors (GMPs) and macrophages (9). As previously discussed, even the network assembling may constitute a useful mechanism to propose novel interactions. This was the case with this model, by envisioning three missing regulations: negative regulation of C/EBP(α) transcription by Foxo1, E2A activation by Ikaros, and Gfi1 positive regulation by Pax-5. Moreover, the model was useful to explore molecular mechanisms of transient induction of the transcription factor C/EBP that down-regulates the transcriptional core of B cell specification and promotes an irreversible trans-differentiation toward macrophages. The theoretical findings complemented the results of a previous experimental report where B cells were transdifferentiated into macrophages by the enforced expression of C/EBP and C/EBP, but without a full understanding on the molecular steps leading to the loss of early and late lymphoid markers and acquisition of myeloid-specific genes (56). Predictions from these models and their perturbations might be useful to unravel the pathobiology of diseases where neoplastic cells concomitantly express myeloid and lymphoid markers (57–61). Also, early branching models may help to deepen the research on plasticity-related processes, such as those suggested to be involved on leukemia lineage switching and relapse (62, 63). It has become of particular interest the integration of microenvironmental cues capable of influencing and regulating transcriptional cores, particularly to approach the two-way feedback between cells and their surrounding microenvironment.

2.3. Microenvironmental Modulation of Lymphoid Differentiation and Plasticity

The applicability of discrete models seems to be simplistic but their scopes are expanding in parallel with the knowledge on cellular heterogeneity and plasticity. Their flexibility for the analysis of biological systems integrated with multiple types of molecular events makes them a useful tool for evaluation of different microenvironments that consider the modulation of genetic and signaling networks. Molecules, such as integrin, cytokine, or antigen receptors, might be included in computational models as they are involved in maintaining particular hematopoietic compartments, enhancing proliferation, regulating apoptosis or migration, or guiding differentiation to either one phenotype or another. As previously mentioned, some of these processes become discrete cellular decisions with a bi-modal behavior as a result of the combined effect of their connectivity in molecular networks and noise (64–66).

Early logical mathematical approaches for modeling lymphocyte behavior upon antigen exposure preceded the development of networks that connected intracellular events regulating the cellular fates of hematopoietic progenitors and lymphocytes (67, 68). However, as the different subtypes of lymphocytes were discovered, efforts focused on the reconstruction of the GRN underlying the emergence of mature phenotypes in response to variable microenvironmental factors under normal and pathological conditions. The first model of lymphoid differentiation branching using a discrete perspective resulted from the transformation of a previous continuous model based on Hill functions describing the polarization of naive Th cells (Th0) into Th1 or Th2 cells (69). The Boolean version proposed by Mendoza (3) integrated 17 nodes and replaced the transcription factor Gata-3 positive self-feedback loop in Yates' model (70) with a more refined functional feedback circuit engaging Gata-3 and interleukin-4 (IL-4) (71). The activation of this functional circuit characterizes the Th2 cell subtype (72, 73). Besides recovering the Th0 polarization into Th1 and Th2, the model was able to describe the transition between Th1 and Th2 attractors by the stimulation with IFN, IL-4, or the combination of IL-12 and IL-18. Later on, the model was extended to include novel transcription factors, cytokines, and signal transduction molecules to describe additional fates to T regulatory (Treg) and Th17 cells (29). More refined molecular data has resulted in the reconstruction of larger versions and their simulations, predicting a larger repertoire of Th cell subsets including Tfh, Th9, Treg, iTreg, Th9, Th17, Th22, and T regulatory Foxp3 independent cells (6, 7, 74).

B and NK cells have been less studied by mathematical modeling. During terminal B cell differentiation in the germinal centers of secondary lymphoid organs, the exposure to particular environmental factors, including the antigen-mediated activation of the B-cell receptor (BCR), defines the transition of the naïve B cell to a memory cell or an antibody-producing plasma cell. This terminal differentiation of B cells has been simulated as a Boolean model that recovered four cellular profiles: naive B cell before and after arriving to the germinal center (GC), memory cell (MC) and plasma cell (PC) (8). The B cell model reproduces not only the expected cellular attractors, but also the transitions with biological significance. Of note, it predicts four interactions that have not been declared experimentally but are suggested through indirect mechanisms: two self-feedback loops involved in Pax5 and Bcl6 activation, the positive regulation of Bcl6 by Pax5, and the inhibition of Pax5 by Irf4.

On the other hand, NK cell biology has been recently approached by a Boolean model providing a CLP attractor that transits toward pro-B, early T progenitor, or three different subtypes of NK attractors, depending on the activation pattern of IL-7, IL-15, and Delta ligand (75). NK cell subsets are characterized by differential expression of the transcription factors T-bet and Eomes. The NK attractor reached after CLP is stimulated with IL-15 activates both transcription factors and correlates with highly cytotoxic NKs both, in humans and mice periphery. On the other hand, perturbation of the CLP attractor with combined activation of IL-15/IL-7 or IL-15/Delta ligand, leads to a T-bet- Eomes⁺ profile correlating with BM NKs or T-bet⁺ Eomes⁻ compatible with liver NKs (76). The incorporation of more transcriptional regulators may lead to new hypothesis about the branching step between NK cells and the more recently described, innate lymphoid cells (ILCs). It has been purposed that CLPs transition to NK lineage may have an intermediate step of a common progenitor for NKs and ILCs with a probable expression of transcription factors shared by both lineages, such as Nfil3 and TOX (77). In contrast to adaptive lymphocytes, the knowledge on transcriptional circuits controlling ILC development remains limited, although their role in the orchestration of immune responses has become of particular interest. ILCs are enriched in mucosal tissues and have been correlated with the progression of allergic, gastrointestinal, and central nervous system inflammatory diseases, like inflammatory bowel disease (IBD) and multiple sclerosis (78, 79). Similar to T lymphocytes, ILCs show plasticity under microenvironmental challenges modifying their cytokine secretion patterns and in consequence, the response exerted by other cells of the immune and adaptive branches (80).

The continuous integration of data, an inevitable process to improve computational modeling of biological systems, leads to the generation of large and complicated networks. To facilitate their analysis, large networks can be subjected to model reduction, a process of iterative removal of particular nodes and redirection of the logical rules that ideally, preserve the reachability of the attractors while keeping the main dynamical properties (29, 81, 82). Model reduction considers that, in a number of cases, a central core of nodes drives the dynamics of other dependent nodes. One of the simplest methodologies to drive model reduction (16) consists in excluding from the steady-state computation those nodes that follow linear downstream pathways. For example, the consecutive rules q₃(t + 1) = q₂(t) and q₂(t + 1) = not q₁(t) may be transformed into q₃ = q₂ = not q₁, so that the state of q₁ automatically determines q₂ and q₃. A more elaborate example would be q₅(t + 1) = q₄(t) and [q₄(t) or q₃(t) or not q₂(t)] which leads to q₅ = q₄ and [q₄ or q₃ or not q₂]; using the Boolean absorption rule a and (a or b) = a, this expression is finally transformed to q₅ = q₄. In this latter case, the steady state of q₅ is merely determined by q₄, independently of the state of q₃ and q₂ which appear in the original rule. Furthermore, model reduction is a useful tool to identify regulatory cores or redundant signal transduction pathways, reduce the states space and obtain qualitative data comparable to experimental results (83–85). An alternative to deal with networks whose size complicates the exhaustive analysis of their state space, consists in the evaluation of cellular transitions assessed with a computational technique known as model checking (6). Model checkers are based on the transformation of states space into graphic or symbolic structures that facilitate verification of properties and trajectories, allowing fate mapping of all possible cell transitions and emerging as a potent predictive tool for cellular plasticity under multiple microenvironmental contexts.

The role of the microenvironment in lymphoid differentiation is successfully implemented in the reviewed models by considering the hypothesis that cytokines are either absent or present, and do not care about graded availability. Other models integrate assumptions to simulate signal processing and propagation using a discrete model, such as the models of the downstream events occurring after the activation of the T-cell receptor (TCR) (11, 12). Saez-Rodriguez et al., based a Boolean model in a large network of 94 nodes and considered that some signaling events occur in a different timescale, so that logical rules were updated in a first and second wave depending on the molecular nature of the event. From the attractors resulting after the simulation, the authors made predictions about the signaling cascade activated by the receptor engagement and confirmed them experimentally. The implementation of two updating waves is a way to recognize that the cellular events occur in different timescales, for example biochemical reactions occurring in the cytoplasm (e.g., molecular inhibition by phosphorylation) are faster than the transcriptional modulations (e.g., transcription factor translocating to the nucleus and binding to a gene promoter that will be activated). Even though their utility, it is necessary to recognize that Boolean models are sometimes insufficient, particularly when there is enough data about the continuous concentration of a biomolecule determinant for the process that is being modeled.

The study of chronic diseases has strongly influenced the understanding of how slight changes derive in the complete perturbation of complex biological systems. If it were desired to simulate the way in which the progressive accumulation of pro-inflammatory factors in the intestinal tract perturb the proportions of T cell populations, the use of Boolean models would be of very limited use to investigate the transitory stages between the healthy attractor and a pathological attractor, like in IBD (78).

3. Modeling of Continuous Variables to Study Lymphocyte Diversity

Modeling lymphoid cells production or activation may require the integration of molecules involved in dosage-dependent effects, as is the case of ligand-receptor affinities, cytokine gradients and even some transcription factors like C/EBP and PU.1 (9, 56). As suggested by the number of publications, continuous mathematical models are the most recurrent tool for the study of lymphocyte development and response and are useful tools to evaluate population dynamics and receptor repertoire (86–89).

However, most time parameters are fitted to experimental data without a deep understanding of molecular mechanisms, unless enough kinetic and biochemical information is available. Some cellular processes involving dosage variations may still be simulated with discrete approaches using multi-valued models or probabilistic Boolean networks, but there exist other alternatives to integrate discrete and continuous molecular events like the construction of hybrid and fuzzy models. On one side, hybrid models have been applied to simulate the activation of Th and B lymphocytes by DCs, and their subsequent departure from the lymph node. The cellular entities and the replication steps were modeled in terms of discrete variables, while the migration was simulated by means of differential equations involving continuous variables and parameters (e.g., chemokines concentration and diffusion, cellular velocity) (90). On the other side, Boolean models may be transformed to continuous systems using fuzzy logic (5, 8, 15–17, 91, 92). These approaches may be useful to use existing GRN of lymphoid differentiation and activation to model complex scenarios that involve intercellular communication among immune cells, interaction of immune cells with normal or pathologic tissue, and immune cell population transitions in response to microenvironment remodeling.

3.1. Dosage Variations in Multi-Valued and Probabilistic Models

The molecular pathways participating in TCR signaling have been successfully modeled with a set of differential equations. The first step for T lymhpocytes activation involves a process known as ligand discrimination that differentiates between weak and strong binding antigens. After TCR engages with peptides processed and expressed on the surface of antigen-presenting cells, a well-regulated discrimination between self and non-self antigens is triggered. The simulation of TCR activation as a continuous model suggested that the MAPK cascade is the responsible for this discriminatory engagement process. A negative feedback loop that modulates the TCR response until an ERK activation threshold is reached may take place, resembling a bimodal behavior (93). The model was expanded to answer the question of how stochastic variations of protein expressions among a clonal population of CD8 T cells could affect their responsiveness. Variations on the expression of CD8 and two components of the MAPK signaling pathway, ERK-1 and SHP-1, generate dispersion in responsiveness among individual cells, but the co-regulation of CD8 and SHP-1 restrain the phenotypic variability (94). It was later discovered that the ligand discrimination process influences T cell differentiation to Treg or Th phenotypes through the downstream modulation of PTEN and Akt/mTOR signaling pathways (95, 96). To represent a weak or strong ligand affinity, a multi-valued model was useful allowing three possible activation levels of TCR and PI3K nodes (off = 0, low = 1, and high = 2). The computational simulations of the model corroborated that low TCR signal favors Treg differentiation, while a stronger signal result in the induction of Th profile (97). Additionally, varying the number of rounds or time-steps for TCR activation, as an approach for ligand binding lifetimes, showed that the Th phenotype is more rapidly stabilized than a Treg profile, suggesting that the transition from naïve to Treg cells is less direct than the Th differentiation. The generation of Treg cells goes through intermediate stages during which the secretion of IL-12 is promoted and activates the PTEN signaling pathway that enables Foxp3 permanent activation (97). Under high TCR signaling, Foxp3 is transiently activated but further turned off by mTOR pathway, while the Akt-dependent regulation of T cell fate choice is also dependent on the differential phosphorylation of additional proteins (98). There are ongoing studies focused on the blockage of TCR signaling by some pathogens like Yersinia pseudotuberculosis (99).

To deepen in the composition of the microenvironmental patterns affecting the diversity of T lymphocytes, a probabilistic Boolean control network (PBCN) was developed for simulation of all possible microenvironments combining nine external signals including TCR activation, TGF-β and IFNγ cytokines, and six interleukines. In contrast with conventional Boolean models, PBCNs contemplate activation probabilities as an approach to input dosages, increasing the range of testable microenvironments (74). Experimental research on T lymphocytes diversity has led to the discovery of intermediate phenotypes that co-express lineage-specific transcription factors from more than one T cell subset, such as Th1-Th2 and Th1-Th17 cells identified on bacterial and parasitic infection (100–102). Through a sensitivity analysis of the PBCN the minimum microenvironment requirements have been identified, on composition and dosage, for the description of each of the 10 T cell profiles. In addition, they have been used to predict the way in which different input patterns influence the internal balance determining the phenotype of canonical and complex cellular profiles, such as cells with mixed phenotypes. With a continuous model constructed to simulate iTreg-Th17 differentiation, Hong and collaborators reported a double expressing phenotype with either regulatory or dual (regulatory and proinflammatory) functions in vivo. This mixed phenotype is suggested to be a stable state reached from the transition of single-expressing cells, iTreg and Th17. Th17 and iTreg cells are able to produce TGF- which may either increase the percentage of both types of cells, or induce the transition from single-expressing to double-expressing cells. The iTreg-Th17 model was also used to analyze how different concentrations of TGF- influence the rate of co-existing cellular subtypes, making evident that priming factors not only drive differentiation events, but also promote cell heterogeneity (103).

The models presented in this section have different limitations. The multi-level and probabilistic models do not allow the integration of temporal hierarchies in the events involved in the biological system of interest, particularly important when modeling more than one type of cellular processes. The continuous model includes a limited number of components depending on the availability of kinetic parameters or enough information to establish assumptions. As an alternative, fuzzy logic can merge large transcriptional regulatory networks participating in cell differentiation and plasticity, with qualitative knowledge about the kinetics of signaling pathways involved in the transduction of microenvironmental variations, for example, events proceeding relatively faster than others, or ligands binding to receptor above other ligands.

3.2. Continuous Simulation of Discrete Differentiation Networks

Extracellular signals and some intracellular components are continuous variables and their adequate representation in mathematical models may determine the simulation of lymphoid cellular fates like differentiation, phenotypic transitions and activation. The transformation of discrete models to a set of differential equations is useful to identify additional attractors and unstable states with biological relevance. In a comparison between Boolean and continuous simulation of a B cell terminal differentiation network, the continuous counterpart provided three additional stable states with intermediate values of Bcl-6 and/or Irf4; however, only one of them was comparable with a previous reported phenotype that may correspond to the centrocytes found during the germinal center selection (8, 104). This intermediate phenotype together with centroblasts, are particularly important in the study of follicular lymphomas characterized by an accumulation of cells unable to reach terminal differentiation stages.

In comparison with Boolean models, the computational simulation of continuous fuzzy models is simpler and in consequence faster, thus allowing the integration of independently developed BRN without caring a lot about the number of resultant equations. An example is the T/B lymphoid differentiation model of 81 equations representing cytokines and transcription factors that lead to ten attractors with Th0, Th1, Th2, Th17, Treg, cytotoxic T lymphocyte, DP T lymphocytes, CD8 T naive, naive B cell, and PC profiles (105). The attractors obtained by the continuous model show a higher compatibility with experimental data than previous discrete models. In this case, all the attractors display intermediate values for Ikaros, Gfi1, and PU1. For each of the three transcription factors there exists strong evidence that associates this intermediate expression with the delimitation toward lymphoid lineage during hematopoietic differentiation (106–108). The intermediate modulation of PU1 and Ikaros was also reproduced with a different continuous model of B-lymphocyte lineage commitment, evidencing their participation in the transcriptional core that reproduces the irreversible transition from LMPP to lineage restricted progenitors expressing IL-7R (109).

An additional application of fuzzy logic models is the simulation of virtual cultures where independent GRNs, representing multiple cells, may interact with a microenvironment expressing graded and dynamic concentrations of cytokines. A virtual culture of T lymphocytes was proposed by Mendoza to evaluate the evolution of 100 cells with an initial Th0 configuration after being stimulated with IFN, I-4, TGF alone, or TGF in combination with IL-6. The phenotype of each cell was determined by the activation state of each of the 36 nodes integrating the internal Th differentiation network, in turn, regulated by 11 cytokines produced depending on each cellular profile (Th0, Treg, Th1, Th2, and Th17). The produced cytokines involved endocrine and paracrine signaling to evaluate the final balance of the T lymphocyte subpopulations arising from different types of stimulus (91). This particular implementation is computationally expensive, but represents a more realistic approach to analyze the interaction between heterogeneous populations of immune cells susceptible to transit among phenotypes, including dynamic secretion patterns that influence the composition of the microenvironment.

3.3. From Discrete to Continuous Using Fuzzy Logic

A more realistic approach must considerate that the expression levels, concentrations, and parameters of biological systems may take any value within a continuous range limited only by functionality constraints. In this case, the discrete dynamic mapping given by Equation (1) may be generalized by introducing a set of ordinary differential equations (ODEs) for the rate of change of the expression level of the network components. For k-th node, this is written as

$\begin{array}{ll}\frac{d{q}_k}{dt}=\mu \left[w_k(q_1,\dots ,q_n)\right]-{\alpha }_k{q}_k.& (2)\end{array}$

Here, μ[w_k] is an input function that expresses a continuous realization of the Boolean rule w_k (see below), while α_k is a decay rate. In this scheme, the equilibrium states of the system are defined by the steady-state condition dq_k/dt = 0, which leads to

$\begin{array}{ll}{q}_k^s=\frac{1}{{\alpha }_k}\mu \left[w_k(q_1^s,\dots ,q_n^s)\right],& (3)\end{array}$

where the superindex s denotes the steady-state value. A straightforward consequence of this is that the expression level of node k is strongly dependent on its decay rate. In the case α_k > 1, a node will be under-expressed with respect to the value attained for α_k = 1; in particular, for α_k ≫ 1, the expression of that node will be completely inhibited: $q_k^s\to 0$. The converse also holds: if α_k < 1, a node will be relatively over-expressed [it must be noticed that a decay rate α_k < 1 may lead to a steady expression value q_k > 1. Although in fuzzy logic the values of the variables are assumed to be constrained to the interval 0 ≤ q_k ≤ 1, values >1 are not excluded by the formalism, and it is a matter of convenience the range in which the variables are defined (110)]. It follows that modifications of the characteristic decay rates of network components may alter the steady expression patterns arising from the nodes interactions. This may be interpreted as a modulation of the metaphorical or Waddington's epigenetic landscape which eventually may lead to transitions between attractors associated to different cell fates. This approach has been formerly employed, for example, to model plastic phenotype changes in T CD4⁺ lymphocytes (92).

The translation of the interactive Boolean rules to the continuous domain may be accomplished by considering an approach based on fuzzy logic. Fuzzy logic is a theory aimed to provide formal foundation to approximate reasoning with applications in physical, biomedical, and behavioral sciences. It is characterized by a graded approach (110–112), so that the degree to which an object exhibits a given property is specified by a membership (or characteristic) function μ[w_k], with truth values ranging from total falsity, μ[w_k] = 0, to totally true, μ[w_k] = 1. For example, the property of “being a good person” implies that there is a set of persons that share certain characteristics with no definite boundary. Fuzzy logic satisfies an axiomatic similar to the implied in Boolean logic, except for the identity principle, meaning that the principle of no-contradiction does not hold. Thus, although seemingly paradoxical, a proposition w and its negation 1 − w may be simultaneously true. For example, the assertion “he was not a good, but not bad guy” has a meaning in language theory. In biological systems, fuzzy propositions may describe cases in which a cell displays an intermediate expression pattern that does not necessarily belong to a specific phenotype. That is the case of individuals with food allergies, in which Treg cells produce IL-4, which is a characteristic usually ascribed to Th2 cells. Similarly, diseases like rheumatoid arthritis or colorectal cancer are associated to the expression of IL-17⁺Foxp3⁺ Treg cells or RORγt + Foxp3⁺ Treg cells, respectively. The absence of no-contradiction is formally expressed by the equation w = 1 − w, with solution w = 1/2. It follows that the value w ≡ w_thr = 1/2 may be interpreted as a threshold between falsity and truth.

Similar to the Boolean approach, in the continuous regime the network regulatory interactions are characterized by fuzzy logic propositions denoted here as w_k[q₁(t), …, q_n(t)]. They are either inferred from experimental observations or suggested by inner consistency requirements. In fact, a translation scheme from the discrete to the continuous scenario may be straightforwardly implemented translation by replacing the Boolean connectives AND, OR, and NOT, for its fuzzy counterparts. In fact, the definition of fuzzy connectives is not unique, and a number of different alternatives not entirely equivalent, have been proposed. In the following table we present Zadeh's original proposal (111) and a probabilistic-like scheme (110):

Both schemes satisfy the modified Boolean axiomatics discussed above. However, the probabilistic-like scheme leads to continuously differentiable expressions if q and p are differentiable. This is a desirable condition when dealing with ODEs systems. Furthermore, it shows the same properties as joint probabilistic distributions for independent variables, so that probabilistic statements may be directly translated into fuzzy propositions.

An example of translation from the Boolean to the fuzzy framework is

$\begin{array}{l}W\left[p,q,r\right]=(q\vee p)\wedge \neg r\to w\left[p,q,r\right]\\ =(q+p-q·p)·(1-r).\end{array}$

Continuous logical propositions can be used to construct an explicit expression of the characteristic function μ[w_k]. In the discrete Boolean approach, this function would be equivalent to a step Θ function:

$\mu [w_k]\to \Theta [w_k-1/2]=\{\begin{array}{cc}0& \text{if\hspace{0.17em}}{w}_k<\text{1/2};\\ 1& \text{if\hspace{0.17em}}{w}_k>\text{1/2}.\end{array}$

In the continuous approach this behavior may be approximated by a characteristic function with a sigmoid structure that gradually changes from a null to a unit value. Many functions share this property. An expression employed in a number of applications of fuzzy logic in systems biology is the logistic function:

$\begin{array}{ll}\mu \left[w_k\right]=\frac{1}{1+exp\left[-\beta (w_k(q_1,\dots ,q_n)-w_{thr})\right]}& (4)\end{array}$

Here, w_thr is a threshold value such that if w_k > w_thr, then w_k tends to be true (or expressed). Usually w_thr = 1/2. The parameter β is a saturation rate that measures the pace of the transit from an unexpressed to an expressed state, as displayed in Figure 3. We observe in the figure that this pace is gradual for small β, and steep for large β. In the case β ≫ 1, μ[w_k] → Θ[w_k − w_thr]. This latter result, together with the steady-state condition given by Equation(2), implies that in the case is another manifestation of the robustness of the qualitative predictions generated by the fuzzy approach. A related result is that in the limit β ≫ 1 and α_k = 1 for every network interaction, then the steady-state condition given by Equation (7), guarantees that the set of fixed-point attractors resulting in the Boolean and fuzzy approaches coincide by construction. On the contrary, the corresponding sets of periodic attractors usually differ.

FIGURE 3

Figure 3. Fuzzy networks. (A) Characteristic distribution $\mu \left[w_k-w^{thr}\right]$ for a threshold value w^thr = 1/2 of the logical proposition w_k, and different values of the saturation parameter β. In the case β ≫ 1, the characteristic distribution becomes a step-like distribution. (B) Fuzzy modeling of the transcriptional core regulating B cell to plasma cell (PC) differentiation using three different saturation values (β = 8, 15, 60), w^thr = 1/2 and the decay rate for each component α = 1. For the initial state of network all nodes were considered inactive, except AP-1, the promoter of the PC differentiation. When β = 60, a full B cell attractor is reached with no final expression of AP-1 or Blimp-1.

It may be argued that the predictions obtained in the fuzzy formalism may depend on the specific form of the characteristic function μ[w_k]. In fact, there are multiple expressions employed for example, in engineering applications and control theory, such as triangular, trapezoidal, or Gaussian functions (113). However, the logistic structure of μ[w_k] considered in this review may be derived, rather than postulated, by introducing the concept of Shannon's information entropy (work in preparation). This is related with the number of independent ways in which a logical proposition may acquire partial values of truth for fixed values of the parameters β and w_thr. In other words, the more general expression involving the least number of assumptions concerning a graded approach to truthiness of a fuzzy proposition is the logistic distribution. Interestingly, the mathematical formalism associated to fuzzy regulatory networks including the description of logical rules with a logistic structure is formally equivalent to that employed in the computation of neural circuits in the theory of neural networks (114).

Another useful (and equivalent) representation of the characteristic function may be derived by considering that the expression levels of biological variables, such as the concentrations or the affinities of a given molecule, may show variations involving several orders of magnitude. In that case, it may be convenient to introduce in the description the logarithms of the corresponding quantities. This is easily performed by means of the change of variable w_k = ln x_k and substituting this into Equation (8), leading to the well-known expression for the Hill function:

$\begin{array}{ll}\tilde{\mu }\left[x_k\right]=\frac{x_k^{\beta }}{x_{thr}^{\beta }+x_k^{\beta }},& (5)\end{array}$

where the parameter x_thr represents the value at which $\tilde{\mu }\left[x_k\right]$ acquires half its maximum value. The Hill function and its negation $1-\tilde{\mu }\left[x_k\right]$ display both a sigmoid shape and have been widely employed in the modeling of biochemical, physiological, and pharmacological processes. A paradigmatic example is the set of non-linear differential equations

$\begin{array}{ll}\frac{d{x}_k}{dt}=\frac{a_k{x}_k^{\beta }}{x_{thr}^{\beta }+x_k^{\beta }}-\frac{b_k{x}_{thr}^{\gamma }}{x_{thr}^{\gamma }+x_k^{\gamma }},& (6)\end{array}$

describing, for example, the simultaneous binding (unbinding)of β (γ) ligands to (from) a single receptor. This latter representation has been employed in the construction of a GRN that characterizes fate decisions and reprogramming signaling pathways of pancreatic cells (115). Although this latter model was not built within the fuzzy logic approach, we observe that in this and numerous instances a formal equivalence may be established by a convenient re-scaling of the variables and parameters involved in the description.

3.4. Self-Organization and Time Ordering

To describe the transitions between distinct steady states, in conjunction with fuzzy logic elements, general concepts of theory of non-equilibrium phase transitions and self-organization are highly relevant to consider. The adaptation of that theory to the fuzzy logic modeling scheme allows a sound description of the transitions between the different disease stages. In the description the transitions between steady states it is important to contemplate that differentiation from a multipotent stem or progenitor state to a mature cell is an essentially irreversible process, and that the associated changes in gene expression patterns exhibit time-directionality. Whereas, in equilibrium systems time-irreversibility is a direct reflection of the second law of thermodynamics, the cell's gene regulatory network represents a system far from thermodynamic equilibrium. These problems have been contemplated by the theory of cooperative phenomena, non-equilibrium phase transition and self-organization (116). Accordingly, cooperative phenomena arise from non-linear interactions of a large number of elementary subsystems (represented here by the fuzzy logic rules), leading to the emergence of organized patterns or phases. The theory relies upon two main concepts, the existence of ordering and control parameters. The order parameters are those variables that mainly drive the system organization, while the control parameters are variables whose value determines which of the possible organizations is actually realized. In the case of thermodynamic systems, an order parameter would be the density, which defines an aggregation state, such as liquid, solid, or gas. These states may somehow “compete,” in the sense that one or other may prevail depending on the value of an external control parameter, such as the temperature of the system, for fixed values of pressure and volume. In the context of fuzzy GRNs, the order parameters are the activation patterns that specify the different cell phenotypes, determined in turn, by the activation state of central nodes or functional moduli of the GRN, while the control parameters are those involved either in the logic rules, or those characterizing the decay rates {α_k}. This latter set is of prime importance. Given that α_k = 1/τ_k, where τ_k is a characteristic expression time, the set {α_k} implicates a hierarchy for the temporal expression of the GRN components. By assuming that an ordering α₁ > α₂ > … can be constructed, this procedure induces an associated ordering τ₁ < τ₂ < …. As in the thermodynamic example, the phenotypic landscape (or state space) may be explored by varying each of the control parameters α_k, while maintaining fixed the rest. As a consequence, transitions between different ordered phases may be induced by modifications of the control parameters. This is similar to the description of chemical reactions in the reaction coordinate space, where the substrate and product states are separated by an activation energy barrier; when an enzyme is added the activation barrier is lowered, and the chemical reaction ensues. In Waddington's landscape context, this mechanism may be interpreted as alterations of the peaks separating the valleys, allowing the exploration of the landscape and transit between valleys. This kind of description has been employed in the modeling of the long-term progression of type-2 diabetes based on a GRN for pancreatic beta cells. In this case, the organization patterns correspond to states identified with health, metabolic syndrome, or manifest diabetes. The alteration of decay rates of key cellular components involved in inflammatory and metabolic pathways lead the transitions between different disease stages.

An important consequence of establishing a time ordering, is that the system dynamics may discriminate among “slow” and “rapid” variables and it may be shown that the main dynamics is driven by “slow” variables, while the “rapid” variables adapt almost instantaneously to the environment defined by the “slow” ones. It turns out that the relaxation times of the order parameters are usually much longer than those of irrelevant variables and thus work like control parameters of the system. Irrelevant variables decay rapidly to a steady state, so that they may be effectively eliminated from the overall dynamics. In this view, the order parameters define the general features of the system, including the final expression patterns associated to a set of initial conditions, while less relevant variables adapt their values to the instructions dictated by the order parameters. This property may be relevant in the study of multifactorial diseases, since it could help in the identification of variables that constitute a target for the development of therapies.

Another element that may be relevant in the study of transitions between steady states is the consideration of extrinsic and intrinsic noise, i.e., the existence of random interactions inherent to every biological system. Depending on its intensity, the existence of noise may drastically alter the predictions yielded by the deterministic formalism considered before, especially at bifurcation points of the landscape, where noise may accelerate a transition rate between neighbor attractors. In the chemical reaction analogy, this is equivalent to adding heat to the process. The action of noise may be incorporated in the fuzzy logic approach by assuming that this is characterized by a stochastic variable ξ(t), with zero mean 〈 ξ(t) 〉 = 0, and statistical dispersion given by 〈 ξ(t) ξ(t′) 〉 = DG(t − t′). Here, D is a diffusion coefficient, and G(t − t′) is a function that characterizes the duration of the self-correlation of the variable ξ. The case in which G is a Dirac delta, i.e., a sharply peaked distribution only for t = t′, and null elsewhere, corresponds to a white noise with no-memory effects.

The fuzzy stochastic dynamics (16) can be described by a Langevin equation (116, 117):

$\begin{array}{ll}\frac{d{q}_k}{dt}=\mu \left[w_k(q_1,\dots ,q_n)\right]-{\alpha }_k{q}_k+{\xi }_k(t),& (7)\end{array}$

with steady states given by the mean value $\langle q_k^s\rangle =\langle \mu \left[w_k^s\right]\rangle /{\alpha }_k$. In the same way as in the deterministic approach, the parameters α_k control the relative heights of peaks and valleys in the mean epigenetic landscape. In the case of small noise (D ≪ 1) the time-dependent solutions are composed by the mean path $\langle q_k^s(t)\rangle$ plus random fluctuations around this path, similarly as dust particles driven by a gentle breeze. The Langevin formalism was implemented by Zhou et al. (115) by means of a GRN addressed to study the processes of differentiation and cell fate reprogramming in pancreatic cells. They show that it is possible to recapitulate the observed attractors of the exocrine and β, δ, α endocrine cells and to predict which gene perturbation can result in a desired lineage reprogramming.

A related approach rests upon a probabilistic or quasi-potential landscape (118, 119). In this case, it is not the ensemble of stochastic trajectories q_k(t), but their probability distribution P[q_k(t)] what constitutes the central concept. One may envisage an epigenetic landscape in which the maximal expression probabilities lie over the deepest (or wider) attraction basins, while the minimal probabilities lie over the hills' tops. Thus, the probabilistic landscape corresponds to an inverted realization of the epigenetic landscape. It can be shown that the probability distribution P[q_k(t)] satisfies the Fokker-Planck (FP) diffusion equation (116, 117), and that the information contained in this formalism is equivalent to that inherent to the Langevin approach. It has been proposed by Wang et al. that the genetic circuitry connections in the quasi-potential landscape imposes the arrow of time in stem cell differentiation, so that the generic asymmetry of barrier heights indicates that the transition from the uncommitted multipotent state to differentiated states is inherently unidirectional.

4. Lymphocyte Involvement in Chronic Diseases: Cellular Diversity and Pathological Feedbacks

The logical framework has also been applied to the simulation of signaling pathways involved in lymphoid related-diseases, like acute lymphoblastic and T cell large granular lymphocyte (T-LGL) leukemia. In the first case, it was predicted that a proinflammatory microenvironment may induce instability in two molecular axes responsible for the retention of hematopoietic progenitor cells within regulatory bone marrow niches (120). In the second case, the model helped to decipher some of the molecular mechanisms that promote survival in T-LGL leukemia cells (121). Both models integrated microenvironmental factors with signaling pathways participating in cellular fate decisions, and in both cases the role of the pro-inflammatory NFkB pathway emerged as important player in the pathogenesis.

Few mathematical models have managed to simulate the dynamic communication between lymphocytes and microenvironment, considering that the feedback loops between both systems are key to modulate immune responses, although the in vivo regulation of both systems is more complex due to influence of neighbor tissues and endocrine signals. The perturbation or inadequate coupling of the regulatory interactions among systems have been suggested to trigger inflammation in multiple chronic diseases. For many years the study of pro-inflammatory conditions was focused on the identification of cytokines as biomarkers or target for adjuvant therapies. With the advances on immunotherapy, the study of immune cells as active participants in chronic diseases development and progression has become of great importance because they represent therapeutic targets with less co-lateral effects than conventional therapies.

Recently, it has been probed that epigenetic landscape approach is useful for the in silico analysis of health to pathogenic progression (122), such as the epithelial-mesenchymal transition and the induction to migratory phenotype induced after chronic pro-inflammatory conditions, offering a tool to delve deeper into transition stages important for early diagnosis (123–126). Computational modeling of epithelial-mesenchymal transition induced by pro-inflammatory cues has suggested an intermediate stage with a senescent profile (125). The process of epithelial malignant transformation is promoted, among other factors, by TGFβ secreted by CD8 and Treg cells, and TNFα produced by macrophages and pro-inflammatory T cells (127, 128). Importantly, CD8 T lymphocytes have been purposed as players in the promotion of aggressive features in breast cancer tumorigenesis (129); but using CD8 T cells as therapeutic targets implicates affecting one of the most important defenses toward infections, so research about the regulatory networks underlying T cell polarization in dynamic feedback with epithelial cells open new opportunities for the development of more precise therapies by simulating multiple or all the possible perturbations in an integral network as also suggested for breast cancer therapy (130).

The same approach is applicable for the study of emergent attractors from many other networks associated to chronic diseases, for example, type 2 diabetes described in terms of beta-pancreatic cell (115) and T lymphocyte interacting networks, based on evidence of the participation of different T subpopulations as inductors of local and systemic inflammation (131). A first approach targeting CD4 T cell plasticity in metabolic diseases showed that hyperinsulinemia, a condition associated with metabolic syndrome and early stages of type 2 diabetes, inhibits the generation of T regulatory Foxp3 cells and stabilizes the Th17 attractor (10). Besides type 2 diabetes, the increase of Th17 subpopulation and decrease of T regulatory cells have been linked with the destruction of beta-pancreatic cells in the pathogenesis of type 1 diabetes, an increased risk of breast cancer recurrence in diabetic patients and increased susceptibility to develop colitis (132–134). The modulation of the Th17 subpopulation as a promising treatment of colitis was predicted by computational simulations of a continuous model. With in silico perturbations of the GRN underlying CD4 T cell differentiation it was predicted that the increase of PPARγ in Th17 cells derives in its transition toward an iTreg profile characterized by the upregulation of Foxp3. The in vivo effect of transplanting PPARγ null Th0 lymphocytes was an increased severity and earlier development of colitis in mice. In contrast, pharmacological activation of PPARγ resulted in the induced shift from Th17 to iTreg phenotype that favored colonic tissue reconstitution (135).

The use of integral models of regulatory networks can be also applied to chronic infections. Existing models of infectious diseases and their interconnection with lymphoid regulatory networks are very limited. Even though, one group has reconstructed a logical network to study the intracellular pathways in CD4 T cells affected by the viral proteins during HIV infection. By considering a model composed by 16 viral and 121 CD4 T cell nodes, they predicted new viral-human molecular interactions and obtained conclusions on the signaling flow affecting cellular fate decisions (136).

All chronic processes mentioned above involve multiple developmental stages where different changes in the microenvironment and the cellular composition take place, depending one on the other through feedback loops. With discrete models we can clearly map the stable stages and the transition between them in the presence or absence of particular nodes, while in conventional continuous models it is quite complicated to include as much as transcription regulators are required to simulate cellular transitions of more than one type of cell. So, the transformation of genetic Boolean models using fuzzy logic, is a promising approach to integrate differentiation networks of lymphoid cells and cells from other tissues to construct more accurate models for the study of chronic diseases, as it becomes important the consideration of temporal evolution and graded changes in molecular compositions. Additionally, less frequent inflammatory cells participating in chronic diseases can be included, like in chronic allergic lung disease, where the progressive accumulation of B cells in the lung promotes Th2 responses by the antigen presentation process (137).

Moreover, intermittent or persistent rapid perturbations in chronic and complex diseases, but not during steady states, provoke small and sometimes, cumulative variations within the cells or their environment, including modifications in cytokines secretion patterns, cellular populations proportions, miRNAs expression, etc., that mostly become visible until there is an abrupt transition of the whole system. Of note, fuzzy logic continuous models permit an easy simulation of such periodical and transient signals that are transduced by cell signaling pathways.

The utility of fuzzy models may apply to a small network composed by some of the main components of the NFkB signaling pathway that behave as a damped oscillator during activation with TNFa (Figure 4). The Boolean simulation of the NFkB network generates two different cyclic attractors when TNFa is activated. However, when simulating the network as a fuzzy logic and varying the parameter of the “slow” reaction corresponding to the genetic transcription of IkBa, damped oscillations as observed in Zambrano et al., are recorded (138). Without introducing any specific kinetic information of receptor affinity, phosphorylation kinetics or translocation velocity, the fuzzy model shows the transition from an initially perturbed system toward a stable state with a controlled or regulated NFkB activation. As suggested by Zambrano et al., these approaches may aid to understand normal cell responses but also their behavior in diseases like cancer, where NFkB activity is usually disregulated and out of control, driving to multiple biological consequences including hyperproliferation, cell survival or migration.

FIGURE 4

Figure 4. Fuzzy models to study signaling pathways activation. (A) NFkB network where IKK is stimulated by microenvironmental TNFα. IKK phosphorylates IkBa unrepressing the dimer p65_RelA to allow its translocation to the nucleus. In the nucleus p65_RelA promotes the transcription of IkBA closing a negative feedback loop of the NFkB pathway. (B) The Boolean simulation of the network generates three attractors, two of them are cyclic attractors with TNFα activation. Here, green = 1, red = 0. (C) Activation value of the node in the NFkB network obtained by fuzzy logic simulation. In this case, β = 3, and α varies depending on the type of biochemical event in which each node is involved. Node tIkBa represents a transcriptional event, with decay rate α = 0.2. (D) Figure taken from (138) showing the nuclear to the cytoplasmic GFP intensity (NCI) of three single GFP-p65 cells stimulated with a constant flow of 10 ng/ml of TNFα.

5. Conclusions

Lymphocytes are active participants of many biological processes involved in homeostasis and can evolve concomitantly to tissues transiting through a pathogenic transformation, due to their responsiveness to a large diversity of biochemical signals and their plasticity. In silico experimentation with regulatory networks has shown the potential to identify the underlying mechanisms of feedback loops that participate in the promotion of disease progression or in the establishment of chronic inflammation. Additionally, the adaptation of existing models for the study of lymphocytes diversity in pathogenic contexts using powerful tools like fuzzy logics represents an approach to visualize the global effect of potential immunotherapeutics.

Author Contributions

JE: analysis of published data, discussion of the topic-related information, drafting, and writing the paper. CV and RP: analysis of published data, discussion of the related information, drafting, writing the paper, critical review of the intellectual content.

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

JE is scholarship holder from CONACyT and IMSS. The authors acknowledge financial support from Instituto de Física, Universidad Nacional Autónoma de México.

References

1. Wang RS, Saadatpour A, Albert R. Boolean modeling in systems biology: an overview of methodology and applications. Phys Biol. (2012) 9:1–14. doi: 10.1088/1478-3975/9/5/055001

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bonzanni N, Garg A, Feenstra KA, Schütte J, Kinston S, Miranda-Saavedra D, et al. Hard-wired heterogeneity in blood stem cells revealed using a dynamic regulatory network model. Bioinformatics (2013) 29:80–8. doi: 10.1093/bioinformatics/btt243

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mendoza L. A network model for the control of the differentiation process in Th cells. Biosystems. (2006) 84:101–14. doi: 10.1016/j.biosystems.2005.10.004

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Naldi A, Thieffry D, Chaouiya C. Decision diagrams for the representation and analysis of logical models of genetic networks. Comput Methods Syst. Biol. (2007) 4695:233–47. doi: 10.1007/978-3-540-75140-3_16

CrossRef Full Text | Google Scholar

5. Martínez-Sosa P, Mendoza L. The regulatory network that controls the differentiation of T lymphocytes. Bio Syst. (2013) 113:96–103. doi: 10.1016/j.biosystems.2013.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Abou-Jaoudé W, Monteiro PT, Naldi A, Grandclaudon M, Soumelis V, Chaouiya C, et al. Model checking to assess T-helper cell plasticity. Front Bioeng Biotechnol. (2015) 2:86. doi: 10.3389/fbioe.2014.00086

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Martinez-Sanchez ME, Mendoza L, Villarreal C, Alvarez-Buylla ER. A minimal regulatory network of extrinsic and intrinsic factors recovers observed patterns of CD4⁺ T cell differentiation and plasticity. PLOS Comput Biol. (2015) 11:e1004324. doi: 10.1371/journal.pcbi.1004324

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Méndez A, Mendoza L. A network model to describe the terminal differentiation of B cells. PLOS Comput. Biol. (2016) 12:e1004696. doi: 10.1371/journal.pcbi.1004696

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Collombet S, van Oevelen C, Sardina Ortega JL, Abou-Jaoudé W, Di Stefano B, Thomas-Chollier M, et al. Logical modeling of lymphoid and myeloid cell specification and transdifferentiation. Proc Natl Acad Sci USA. (2017) 114:5792–9. doi: 10.1073/pnas.1610622114

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Martinez-Sanchez ME, Hiriart M, Alvarez-Buylla ER. The CD4⁺ T cell regulatory network mediates inflammatory responses during acute hyperinsulinemia: a simulation study. BMC Syst Biol. (2017) 11:64. doi: 10.1186/s12918-017-0436-y

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Saez-Rodriguez J, Simeoni L, Lindquist JA, Hemenway R, Bommhardt U, Arndt B, et al. A logical model provides insights into T cell receptor signaling. PLoS Comput Biol. (2007) 3:e163. doi: 10.1371/journal.pcbi.0030163

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles ED. A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics (2006) 7:56. doi: 10.1186/1471-2105-7-56

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hunter SD, Peters LE, Wotherspoon JS, Crowe SM. Lymphocyte subset analysis by boolean algebra: a phenotypic approach using a cocktail of 5 antibodies and 3 color immunofluorescence. Cytometry. (1994) 15:258–66. doi: 10.1002/cyto.990150311

PubMed Abstract | CrossRef Full Text | Google Scholar

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

15. Di Cara A, Garg A, De Micheli G, Xenarios I, Mendoza L. Dynamic simulation of regulatory networks using SQUAD. BMC Bioinformatics. (2007) 8:462. doi: 10.1186/1471-2105-8-462

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Villarreal C, Padilla-Longoria P, Alvarez-Buylla ER. General theory of genotype to phenotype mapping: derivation of epigenetic landscapes from N-node complex gene regulatory networks. Phys Rev Lett. (2012) 109:118102. doi: 10.1103/PhysRevLett.109.118102

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Morris MK, Saez-Rodriguez J, Clarke DC, Sorger PK, Lauffenburger DA. Training signaling pathway maps to biochemical data with constrained fuzzy logic: quantitative analysis of liver cell responses to inflammatory stimuli. PLoS Comput Biol. (2011) 7:e1001099. doi: 10.1371/journal.pcbi.1001099

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Waddington CH. The Strategy of the Genes. London: George Allen & Unwin (1957).

Google Scholar

19. Kauffman SAa. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. (1969) 22:437–67. doi: 10.1016/0022-5193(69)90015-0

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Garcia HG, Tikhonov M, Lin A, Gregor T. Quantitative imaging of transcription in living Drosophila embryos links polymerase activity to patterning. Curr Biol. (2013) 23:2140–5. doi: 10.1016/j.cub.2013.08.054

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Aldana M. Boolean dynamics of networks with scale-free topology. Phys D Nonlinear Phen. (2003) 185:45–66. doi: 10.1016/S0167-2789(03)00174-X

CrossRef Full Text | Google Scholar

22. Khanin R, Wit E. How scale-free are biological networks. J Comput Biol. (2006) 13:810–8. doi: 10.1089/cmb.2006.13.810

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Serra R, Villani M, Agostini L. On the dynamics of scale-free Boolean networks. In: Apolloni B, Marinaro M, Tagliaferri R, editors. Neural Nets. Berlin; Heidelberg: Springer (2003). p. 43–9.

Google Scholar

24. Albert R. Scale-free networks in cell biology. J Cell Sci. (2005) 118:4947–57. doi: 10.1242/jcs.02714

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Thomas R. On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. In: Della Dora J, Demongeot J, Lacolle B, editors. Numerical Methods in the Study of Critical Phenomena. Springer Series in Synergetics, Vol. 9. Berlin; Heidelberg: Springer (1981). p. 180–93.

Google Scholar

26. Fauré A, Naldi A, Chaouiya C, Thieffry D, Faure A, Naldi A, et al. Dynamical analysis of a generic Boolean model for the control of the mammalian cell cycle. Bioinformatics. (2006) 22:e124–31. doi: 10.1093/bioinformatics/btl210

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Davidich MI, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PLoS ONE. (2008) 3:e1672. doi: 10.1371/journal.pone.0001672

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ortiz-Gutiérrez E, García-Cruz K, Azpeitia E, Castillo A, Sánchez MdlP, Álvarez-Buylla ER. A dynamic gene regulatory network model that recovers the cyclic behavior of Arabidopsis thaliana cell cycle. PLOS Comput Biol. (2015) 11:e1004486. doi: 10.1371/journal.pcbi.1004486

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Naldi A, Carneiro J, Chaouiya C, Thieffry D. Diversity and plasticity of Th cell types predicted from regulatory network modelling. PLoS Comput Biol. (2010) 6:e1000912. doi: 10.1371/journal.pcbi.1000912

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Bhattacharya S, Conolly RB, Kaminski NE, Thomas RS, Andersen ME, Zhang Q. A bistable switch underlying B-cell differentiation and its disruption by the environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin. Toxicol Sci. (2010) 115:51–65. doi: 10.1093/toxsci/kfq035

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Montecino-Rodriguez E, Leathers H, Dorshkind K. Bipotential B-macrophage progenitors are present in adult bone marrow. Nat Immunol. (2001) 2:83–8. doi: 10.1038/83210

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Iwasaki H, Akashi K. Hematopoietic developmental pathways: on cellular basis. Oncogene. (2007) 26:6687–96. doi: 10.1038/sj.onc.1210754

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hao QL, Zhu J, Price MA, Payne KJ, Barsky LW, Crooks GM. Identification of a novel, human multilymphoid progenitor in cord blood. Blood. (2001) 97:3683–90. doi: 10.1182/blood.V97.12.3683

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ghaedi M, Steer C, Martinez-Gonzalez I, Halim TF, Abraham N, Takei F. Common-lymphoid-progenitor-independent pathways of innate and T lymphocyte development. Cell Rep. (2016) 15:471–80. doi: 10.1016/j.celrep.2016.03.039

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Arinobu Y, Mizuno Si, Chong Y, Shigematsu H, Iino T, Iwasaki H, et al. Reciprocal activation of GATA-1 and PU.1 marks initial specification of hematopoietic stem cells into myeloerythroid and myelolymphoid lineages. Cell Stem Cell. (2007) 1:416–27. doi: 10.1016/j.stem.2007.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Grzywacz B, Kataria N, Kataria N, Blazar BR, Miller JS, Verneris MR. Natural killer-cell differentiation by myeloid progenitors. Blood. (2011) 117:3548–58. doi: 10.1182/blood-2010-04-281394

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Akashi K, Traver D, Miyamoto T, Weissman IL. A clonogenic common myeloid progenitor that gives rise to all myeloid lineages. Nature. (2000) 404:193–7. doi: 10.1038/35004599

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Christensen JL, Weissman IL. Flk-2 is a marker in hematopoietic stem cell differentiation: a simple method to isolate long-term stem cells. Proc Natl Acad Sci USA. (2001) 98:14541–6. doi: 10.1073/pnas.261562798

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Adolfsson J, Månsson R, Buza-Vidas N, Hultquist A, Liuba K, Jensen CT, et al. Identification of Flt3⁺ lympho-myeloid stem cells lacking erythro-megakaryocytic potential. Cell. (2005) 121:295–306. doi: 10.1016/j.cell.2005.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Galy A, Travis M, Cen D, Chen B. Human T, B, natural killer, and dendritic cells arise from a common bone marrow progenitor cell subset. Immunity. (1995) 3:459–73. doi: 10.1016/1074-7613(95)90175-2

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Fathman JW, Bhattacharya D, Inlay MA, Seita J, Karsunky H, Weissman IL. Identification of the earliest natural killer cell-committed progenitor in murine bone marrow. Blood. (2011) 118:5439–47. doi: 10.1182/blood-2011-04-348912

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Karsunky H, Inlay MA, Serwold T, Bhattacharya D, Weissman IL. Flk2⁺ common lymphoid progenitors possess equivalent differentiation potential for the B and T lineages. Blood. (2008) 111:5562–70. doi: 10.1182/blood-2007-11-126219

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Serwold T, Richie Ehrlich LI, Weissman IL. Reductive isolation from bone marrow and blood implicates common lymphoid progenitors as the major source of thymopoiesis. Blood. (2009) 113:807–15. doi: 10.1182/blood-2008-08-173682

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Inlay MA, Bhattacharya D, Sahoo D, Serwold T, Seita J, Karsunky H, et al. Ly6d marks the earliest stage of B-cell specification and identifies the branchpoint between B-cell and T-cell development. Genes Dev. (2009) 23:2376–81. doi: 10.1101/gad.1836009

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Kondo M, Scherer DC, Miyamoto T, King AG, Akashi K, Sugamura K, et al. Cell-fate conversion of lymphoid-committed progenitors by instructive actions of cytokines. Nature. (2000) 407:383–6. doi: 10.1038/35030112

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Constantinides MG, McDonald BD, Verhoef PA, Bendelac A. A committed precursor to innate lymphoid cells. Nature. (2014) 508:397–401. doi: 10.1038/nature13047

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Welner RS, Pelayo R, Nagai Y, Garrett KP, Wuest TR, Carr DJ, et al. Lymphoid precursors are directed to produce dendritic cells as a result of TLR9 ligation during herpes infection. Blood. (2008) 112:3753–61. doi: 10.1182/blood-2008-04-151506

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Welner RS, Pelayo R, Garrett KP, Chen X, Perry SS, Sun Xh, et al. Interferon-producing killer dendritic cells (IKDCs) arise via a unique differentiation pathway from primitive c-kit Hi CD62L⁺ lymphoid progenitors. Blood. (2007) 109:4825–31. doi: 10.1182/blood-2006-08-043810

CrossRef Full Text | Google Scholar

49. Pelayo R, Hirose J, Huang J, Garrett KP, Delogu A, Busslinger M, et al. Derivation of 2 categories of plasmacytoid dendritic cells in murine bone marrow. Blood. (2005) 105:4407–15. doi: 10.1182/blood-2004-07-2529

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Welner RS, Pelayo R, Kincade PW. Evolving views on the genealogy of B cells. Nat Rev Immunol. (2008) 8:95–106. doi: 10.1038/nri2234

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Zlotoff DA, Bhandoola A. Hematopoietic progenitor migration to the adult thymus. Ann N Y Acad Sci. (2011) 1217:122–38. doi: 10.1111/j.1749-6632.2010.05881.x

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Herzog S, Reth M, Jumaa H. Regulation of B-cell proliferation and differentiation by pre-B-cell receptor signalling. Nat Rev Immunol. (2009) 9:195–205. doi: 10.1038/nri2491

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Guo H, Barberi T, Suresh R, Friedman AD. Progression from the common lymphoid progenitor to B/myeloid PreproB and ProB precursors during B lymphopoiesis requires C/EBPα. J Immunol. (2018) 201:1692–704. doi: 10.4049/jimmunol.1800244

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ramos CA, Bowman TA, Boles NC, Merchant AA, Zheng Y, Parra I, et al. Evidence for diversity in transcriptional profiles of single hematopoietic stem cells. PLoS Genet. (2006) 2:e159. doi: 10.1371/journal.pgen.0020159

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Chambers SM, Boles NC, Lin KYK, Tierney MP, Bowman TV, Bradfute SB, et al. Hematopoietic fingerprints: an expression database of stem cells and their progeny. Cell Stem Cell. (2007) 1:578–91. doi: 10.1016/j.stem.2007.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Xie H, Ye M, Feng R, Graf T. Stepwise reprogramming of B cells into macrophages. Cell. (2004) 117:663–76. doi: 10.1016/S0092-8674(04)00419-2

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Scamurra DO, Davey FR, Nelson DA, Kurec AS, Goldberg J. Acute leukemia presenting with myeloid and lymphoid cell markers. Ann Clin Lab Sci. (1983) 13:496–502.

PubMed Abstract | Google Scholar

58. Meyer PR, Krugliak L, Neely S, Levine A, Parker JW, Kaplan B, et al. Acute leukemias with both myeloid and lymphoid surface markers: cytoplasmic alpha-1-anti-chymotrypsin and alpha-1-anti-trypsin as possible indicators of early granulocytic differentiation. Am J Clin Pathol. (1986) 86:461–8. doi: 10.1093/ajcp/86.4.461

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Mirro J, Kitchingman GR, Williams DL, Murphy SB, Zipf TF, Stass SA. Mixed lineage leukemia: the implications for hematopoietic differentiation. Blood. (1986) 68:597–9.

PubMed Abstract | Google Scholar

60. Quesada AE, Hu Z, Routbort MJ, Patel KP, Luthra R, Loghavi S, et al. Mixed phenotype acute leukemia contains heterogeneous genetic mutations by next-generation sequencing. Oncotarget. (2018) 9:8441–9. doi: 10.18632/oncotarget.23878

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Hammond WA, Advani P, Ketterling RP, Van Dyke D, Foran JM, Jiang L. Biphenotypic acute leukemia versus myeloid antigen-positive ALL: clinical relevance of WHO criteria for mixed phenotype acute leukemia. Case Rep Hematol. (2018) 2018:7456378. doi: 10.1155/2018/7456378

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Dorantes-Acosta E, Pelayo R. Lineage switching in acute leukemias: a consequence of stem cell plasticity? Bone Marrow Res. (2012) 2012:406796. doi: 10.1155/2012/406796

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Jacoby E, Nguyen SM, Fountaine TJ, Welp K, Gryder B, Qin H, et al. CD19 CAR immune pressure induces B-precursor acute lymphoblastic leukaemia lineage switch exposing inherent leukaemic plasticity. Nat Commun. (2016) 7:12320. doi: 10.1038/ncomms12320

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Huang CY, Ferrell JE. Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc Natl Acad Sci USA. (1996) 93:10078–83. doi: 10.1073/pnas.93.19.10078

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Kim K, Sauro HM. In search of noise-induced bimodality. BMC Biol. (2012) 10:89. doi: 10.1186/1741-7007-10-89

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Das J, Ho M, Zikherman J, Govern C, Yang M, Weiss A, et al. Digital signaling and hysteresis characterize Ras activation in lymphoid cells. Cell. (2009) 136:337–51. doi: 10.1016/j.cell.2008.11.051

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Kaufman M, Urbain J, Thomas R. Towards a logical analysis of the immune response. J Theor Biol. (1985) 114:527–61. doi: 10.1016/S0022-5193(85)80042-4

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kaufman M, Thomas R. Model analysis of the bases of multistationarity in the humoral immune response. J Theor Biol. (1987) 129:141–62. doi: 10.1016/S0022-5193(87)80009-7

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Yates A, Callard R, Stark J. Combining cytokine signalling with T-bet and GATA-3 regulation in Th1 and Th2 differentiation: a model for cellular decision-making. J Theor Biol. (2004) 231:181–96. doi: 10.1016/j.jtbi.2004.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Ouyang W, Löhning M, Gao Z, Assenmacher M, Ranganath S, Radbruch A, et al. Stat6-independent GATA-3 autoactivation directs IL-4-independent Th2 development and commitment. Immunity. (2000) 12:27–37. doi: 10.1016/S1074-7613(00)80156-9

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Zheng W, Flavell RA. The transcription factor GATA-3 is necessary and sufficient for Th2 cytokine gene expression in CD4 T cells. Cell. (1997) 89:587–96. doi: 10.1016/S0092-8674(00)80240-8

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Kim JI, Ho IC, Grusby MJ, Glimcher LH. The transcription factor c-Maf controls the production of interleukin-4 but not other Th2 cytokines. Immunity. (1999) 10:745–51. doi: 10.1016/S1074-7613(00)80073-4

CrossRef Full Text | Google Scholar

73. Ho IC, Hodge MR, Rooney JW, Glimcher LH. The proto-oncogene c-maf is responsible for tissue-specific expression of interleukin-4. Cell. (1996) 85:973–83. doi: 10.1016/S0092-8674(00)81299-4

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Puniya BL, Todd RG, Mohammed A, Brown DM, Barberis M, Helikar T. A mechanistic computational model reveals that plasticity of CD4⁺ T cell differentiation is a function of cytokine composition and dosage. Front Physiol. (2018) 9:878. doi: 10.3389/fphys.2018.00878

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Liquitaya-Montiel AJ, Mendoza L. Dynamical analysis of the regulatory network controlling natural killer cells differentiation. Front Physiol. (2018) 9:1029. doi: 10.3389/fphys.2018.01029

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Simonetta F, Pradier A, Roosnek E. T-bet and eomesodermin in NK cell development, maturation, and function. Front Immunol. (2016) 7:241. doi: 10.3389/fimmu.2016.00241

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Zhong C, Zhu J. Transcriptional regulatory network for the development of innate lymphoid cells. Mediat Inflamm. (2015) 2015:264502. doi: 10.1155/2015/264502

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Geremia A, Arancibia-Cárcamo CV. Innate lymphoid cells in intestinal inflammation. Front Immunol. (2017) 8:1296. doi: 10.3389/fimmu.2017.01296

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Brown MA, Weinberg RB. Mast cells and innate lymphoid cells: underappreciated players in CNS autoimmune demyelinating disease. Front Immunol. (2018) 9:514. doi: 10.3389/fimmu.2018.00514

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Lim AI, Verrier T, Vosshenrich CA, Di Santo JP. Developmental options and functional plasticity of innate lymphoid cells. Curr Opin Immunol. (2017) 44:61–8. doi: 10.1016/j.coi.2017.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Naldi A, Remy E, Thieffry D, Chaouiya C. Dynamically consistent reduction of logical regulatory graphs. Theor Comput Sci. (2011) 412:2207–18. doi: 10.1016/j.tcs.2010.10.021

CrossRef Full Text | Google Scholar

82. Saadatpour A, Albert R, Reluga TC. A reduction method for Boolean network models proven to conserve attractors. SIAM J Appl Dyn Syst. (2013) 12:1997–2011. doi: 10.1137/13090537X

CrossRef Full Text | Google Scholar

83. Calzone L, Tournier L, Fourquet S, Thieffry D, Zhivotovsky B, Barillot E, et al. Mathematical modelling of cell-fate decision in response to death receptor engagement. PLoS ONE. (2010) 6:e1000702. doi: 10.1371/journal.pcbi.1000702

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Bérenguier D, Chaouiya C, Monteiro PT, Naldi A, Remy E, Thieffry D, et al. Dynamical modeling and analysis of large cellular regulatory networks. Chaos. (2013) 23:025114. doi: 10.1063/1.4809783

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Mendes ND, Lang F, Le Cornec YS, Mateescu R, Batt G, Chaouiya C. Composition and abstraction of logical regulatory modules: application to multicellular systems. Bioinformatics. (2013) 29:749–57. doi: 10.1093/bioinformatics/btt033

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Mehr R. Asynchronous differentiation models explain bone marrow labeling kinetics and predict reflux between the pre- and immature B cell pools. Int Immunol. (2003) 15:301–12. doi: 10.1093/intimm/dxg025

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Shahaf G, Allman D, Cancro MP, Mehr R. Screening of alternative models for transitional B cell maturation. Int Immunol. (2004) 16:1081–90. doi: 10.1093/intimm/dxh109

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Salmon-Divon M, Höglund P, Mehr R. Generation of the natural killer cell repertoire: the sequential vs. the two-step selection model. Bull Math Biol. (2003) 65:199–218. doi: 10.1016/S0092-8240(02)00093-9

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Salmon-Divon M, Höglund P, Johansson MH, Johansson S, Mehr R. Computational modeling of human natural killer cell development suggests a selection process regulating coexpression of KIR with CD94/NKG2A. Mol Immunol. (2005) 42:397–403. doi: 10.1016/j.molimm.2004.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

90. Baldazzi V, Paci P, Bernaschi M, Castiglione F. Modeling lymphocyte homing and encounters in lymph nodes. BMC Bioinformatics. (2009) 10:387. doi: 10.1186/1471-2105-10-387

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Mendoza L. A virtual culture of CD4⁺ T lymphocytes. Bull Math Biol. (2013) 75:1012–29. doi: 10.1007/s11538-013-9814-9

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Martinez-Sanchez ME, Huerta L, Alvarez-Buylla ER, Villarreal Luján C. Role of cytokine combinations on CD4⁺ T cell differentiation, partial polarization, and plasticity: continuous network modeling approach. Front Physiol. (2018) 9:877. doi: 10.3389/fphys.2018.00877

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Altan-Bonnet G, Germain RN. Modeling T cell antigen discrimination based on feedback control of digital ERK responses. PLoS Biol. (2005) 3:e356. doi: 10.1371/journal.pbio.0030356

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Feinerman O, Veiga J, Dorfman JR, Germain RN, Altan-Bonnet G. Variability and robustness in T cell activation from regulated heterogeneity in protein levels. Science. (2008) 321:1081–4. doi: 10.1126/science.1158013

PubMed Abstract | CrossRef Full Text | Google Scholar

95. Sauer S, Bruno L, Hertweck A, Finlay D, Leleu M, Spivakov M, et al. T cell receptor signaling controls Foxp3 expression via PI3K, Akt, and mTOR. Proc Natl Acad Sci USA. (2008) 105:7797–802. doi: 10.1073/pnas.0800928105

PubMed Abstract | CrossRef Full Text | Google Scholar

96. Haxhinasto S, Mathis D, Benoist C. The AKT-mTOR axis regulates de novo differentiation of CD4⁺Foxp3⁺ cells. J Exp Med. (2008) 205:565–74. doi: 10.1084/jem.20071477

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Miskov-Zivanov N, Turner MS, Kane LP, Morel PA, Faeder JR. The duration of T cell stimulation is a critical determinant of cell fate and plasticity. Sci Signal. (2013) 6:ra97. doi: 10.1126/scisignal.2004217

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Hawse WF, Boggess WC, Morel PA. TCR signal strength regulates Akt substrate specificity to induce alternate murine Th and T regulatory cell differentiation programs. J Immunol. (2017) 199:589–97. doi: 10.4049/jimmunol.1700369

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Pasztoi M, Bonifacius A, Pezoldt J, Kulkarni D, Niemz J, Yang J, et al. Yersinia pseudotuberculosis supports Th17 differentiation and limits de novo regulatory T cell induction by directly interfering with T cell receptor signaling. Cell Mol Life Sci. (2017) 74:2839–50. doi: 10.1007/s00018-017-2516-y

PubMed Abstract | CrossRef Full Text | Google Scholar

100. Kullberg MC, Jankovic D, Feng CG, Hue S, Gorelick PL, McKenzie BS, et al. IL-23 plays a key role in Helicobacter hepaticus-induced T cell-dependent colitis. J Exp Med. (2006) 203:2485–94. doi: 10.1084/jem.20061082

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Morrison PJ, Bending D, Fouser LA, Wright JF, Stockinger B, Cooke A, et al. Th17-cell plasticity in Helicobacter hepaticus-induced intestinal inflammation. Mucosal Immunol. (2013) 6:1143–56. doi: 10.1038/mi.2013.11

PubMed Abstract | CrossRef Full Text | Google Scholar

102. Bock CN, Babu S, Breloer M, Rajamanickam A, Boothra Y, Brunn ML, et al. Th2/1 hybrid cells occurring in murine and human strongyloidiasis share effector functions of Th1 cells. Front Cell Infect Microbiol. (2017) 7:261. doi: 10.3389/fcimb.2017.00261

PubMed Abstract | CrossRef Full Text | Google Scholar

103. Hong T, Xing J, Li L, Tyson JJ. A mathematical model for the reciprocal differentiation of T helper 17 cells and induced regulatory T cells. PLoS Comput Biol. (2011) 7:e1002122. doi: 10.1371/journal.pcbi.1002122

PubMed Abstract | CrossRef Full Text | Google Scholar

104. Cattoretti G, Shaknovich R, Smith PM, Jäck HM, Murty VV, Alobeid B. Stages of germinal center transit are defined by B cell transcription factor coexpression and relative abundance. J Immunol. (2006) 177:6930–9. doi: 10.4049/jimmunol.177.10.6930

PubMed Abstract | CrossRef Full Text | Google Scholar

105. Mendoza L, Méndez A. A dynamical model of the regulatory network controlling lymphopoiesis. Biosystems. (2015) 137:26–33. doi: 10.1016/j.biosystems.2015.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

106. Spooner CJ, Cheng JX, Pujadas E, Laslo P, Singh H. A recurrent network involving the transcription factors PU.1 and Gfi1 orchestrates innate and adaptive immune cell fates. Immunity. (2009) 31:576–86. doi: 10.1016/j.immuni.2009.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

107. Mak KS, Funnell APW, Pearson RCM, Crossley M. PU.1 and haematopoietic cell fate: dosage matters. Int J Cell Biol. (2011) 2011:808524. doi: 10.1155/2011/808524

PubMed Abstract | CrossRef Full Text | Google Scholar

108. Zarnegar MA, Rothenberg EV. Ikaros represses and activates PU.1 cell-type-specifically through the multifunctional Sfpi1 URE and a myeloid specific enhancer. Oncogene. (2012) 31:4647–54. doi: 10.1038/onc.2011.597

PubMed Abstract | CrossRef Full Text | Google Scholar

109. Salerno L, Cosentino C, Morrone G, Amato F. Computational modeling of a transcriptional switch underlying B-lymphocyte lineage commitment of hematopoietic multipotent. PLoS ONE. (2015) 10:e132208. doi: 10.1371/journal.pone.0132208

PubMed Abstract | CrossRef Full Text | Google Scholar

110. Novak V, Perfilieva I, Mockor J. Mathematical Principles of Fuzzy Logic. New York, NY: Kluwer Academic (1999).

Google Scholar

111. Zadeh LA. Fuzzy logic and approximate reasoning. Synthese. (1975) 30:407–28. doi: 10.1007/BF00485052

CrossRef Full Text | Google Scholar

112. Dubois D, Prade H. The three semantics of fuzzy sets. Fuzzy Sets Syst. (1997) 90:141–50. doi: 10.1016/S0165-0114(97)00080-8

CrossRef Full Text | Google Scholar

113. Ross TJ. Fuzzy Logic With Engineering Applications. West Sussex: John Wiley (2004).

Google Scholar

114. Hopfield JJ. Neural networks and physical systems with emergent collective computational abilities. Proc Natl Acad Sci USA. (1982) 79:2554–8. doi: 10.1073/pnas.79.8.2554

PubMed Abstract | CrossRef Full Text | Google Scholar

115. Zhou JX, Brusch L, Huang S. Predicting pancreas cell fate decisions and reprogramming with a hierarchical multi-attractor model. PLoS ONE. (2011) 6:e14752. doi: 10.1371/journal.pone.0014752

PubMed Abstract | CrossRef Full Text | Google Scholar

116. Haken H. Synergetics: An Introduction: Nonequilibrium Phase Transitions and Self-organization in Physics, Chemistry, and Biology. Springer (1983).

Google Scholar

117. Haken H. Cooperative phenomena in systems far from thermal equilibrium and in nonphysical systems. Rev Modern Phys. (1975) 47:67–121. doi: 10.1103/RevModPhys.47.67

CrossRef Full Text | Google Scholar

118. Wang J, Xu L, Wang E, Huang S. The potential landscape of genetic circuits imposes the arrow of time in stem cell differentiation. Biophys J. (2010) 99:29–39. doi: 10.1016/j.bpj.2010.03.058

PubMed Abstract | CrossRef Full Text | Google Scholar

119. Wang J, Xu L, Wang E. Potential landscape and flux framework of nonequilibrium networks: robustness, dissipation, and coherence of biochemical oscillations. Proc Natl Acad Sci USA. (2008) 105:12271–6. doi: 10.1073/pnas.0800579105

PubMed Abstract | CrossRef Full Text | Google Scholar

120. Enciso J, Mayani H, Mendoza L, Pelayo R. Modeling the pro-inflammatory tumor microenvironment in acute lymphoblastic leukemia predicts a breakdown of hematopoietic-mesenchymal communication networks. Front Physiol. (2016) 7:349. doi: 10.3389/fphys.2016.00349

PubMed Abstract | CrossRef Full Text | Google Scholar

121. Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, et al. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci USA. (2008) 105:16308–13. doi: 10.1073/pnas.0806447105

PubMed Abstract | CrossRef Full Text | Google Scholar

122. Huang S, Ernberg I, Kauffman S. Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective. Semin Cell Dev Biol. (2009) 20:869–76. doi: 10.1016/j.semcdb.2009.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

123. Cohen DPA, Martignetti L, Robine S, Barillot E, Zinovyev A, Calzone L. Mathematical modelling of molecular pathways enabling tumour cell invasion and migration. PLoS Comput Biol. (2015) 11:e1004571. doi: 10.1371/journal.pcbi.1004571

PubMed Abstract | CrossRef Full Text | Google Scholar

124. Montagud A, Traynard P, Martignetti L, Bonnet E, Barillot E, Zinovyev A, et al. Conceptual and computational framework for logical modelling of biological networks deregulated in diseases. Brief Bioinform. (2017). doi: 10.1093/bib/bbx163. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

125. Méndez-López LF, Davila-Velderrain J, Domínguez-Hüttinger E, Enríquez-Olguín C, Martínez-García JC, Alvarez-Buylla ER. Gene regulatory network underlying the immortalization of epithelial cells. BMC Syst Biol. (2017) 11:24. doi: 10.1186/s12918-017-0393-5

PubMed Abstract | CrossRef Full Text | Google Scholar

126. Fumiã HF, Martins ML. Boolean network model for cancer pathways: predicting carcinogenesis and targeted therapy outcomes. PLoS ONE. (2013) 8:e69008. doi: 10.1371/journal.pone.0069008

PubMed Abstract | CrossRef Full Text | Google Scholar

127. Bates RC, Mercurio AM. Tumor necrosis factor-α stimulates the epithelial-to-mesenchymal transition of human colonic organoids. Mol Biol Cell. (2003) 14:1790–800. doi: 10.1091/mbc.e02-09-0583

PubMed Abstract | CrossRef Full Text | Google Scholar

128. Mahic M, Henjum K, Yaqub S, Bjørnbeth BA, Torgersen KM, Taskén K, et al. Generation of highly suppressive adaptive CD8⁺CD25⁺FOXP3⁺ regulatory T cells by continuous antigen stimulation. Eur J Immunol. (2008) 38:640–6. doi: 10.1002/eji.200737529

PubMed Abstract | CrossRef Full Text | Google Scholar

129. Santisteban M, Reiman JM, Asiedu MK, Behrens MD, Nassar A, Kalli KR, et al. Immune-induced epithelial to mesenchymal transition in vivo generates breast cancer stem cells. Cancer Res. (2009) 69:2887–95. doi: 10.1158/0008-5472.CAN-08-3343

PubMed Abstract | CrossRef Full Text | Google Scholar

130. Gómez Tejeda Zañudo J, Scaltriti M, Albert R. A network modeling approach to elucidate drug resistance mechanisms and predict combinatorial drug treatments in breast cancer. Cancer Conv. (2017) 1:5. doi: 10.1186/s41236-017-0007-6

PubMed Abstract | CrossRef Full Text | Google Scholar

131. Xia C, Rao X, Zhong J. Role of T lymphocytes in type 2 diabetes and diabetes-associated inflammation. J Diabet Res. (2017) 2017:1–6. doi: 10.1155/2017/6494795

PubMed Abstract | CrossRef Full Text | Google Scholar

132. Emamaullee JA, Davis J, Merani S, Toso C, Elliott JF, Thiesen A, et al. Inhibition of Th17 cells regulates autoimmune diabetes in NOD mice. Diabetes. (2009) 58:1302–11. doi: 10.2337/db08-1113

PubMed Abstract | CrossRef Full Text | Google Scholar

133. Nicholas D, Strissel KJ, Andrieu G, Tran AH, Nikolajczyk BS, Denis GV. The type 2 diabetes-associated Th17 cytokines IL-21 and IL-22 promote breast cancer cell survival. J Immunol. (2016) 196:124.56.

Google Scholar

134. Xin L, Jiang TT, Chaturvedi V, Kinder JM, Ertelt JM, Rowe JH, et al. Commensal microbes drive intestinal inflammation by IL-17-producing CD4⁺ T cells through ICOSL and OX40L costimulation in the absence of B7-1 and B7-2. Proc Natl Acad Sci USA. (2014) 111:10672–7. doi: 10.1073/pnas.1402336111

PubMed Abstract | CrossRef Full Text | Google Scholar

135. Carbo A, Hontecillas R, Kronsteiner B, Viladomiu M, Pedragosa M, Lu P, et al. Systems modeling of molecular mechanisms controlling cytokine-driven CD4⁺ T cell differentiation and phenotype plasticity. PLoS Comput Biol. (2013) 9:e1003027. doi: 10.1371/journal.pcbi.1003027

PubMed Abstract | CrossRef Full Text | Google Scholar

136. Oyeyemi OJ, Davies O, Robertson DL, Schwartz JM. A logical model of HIV-1 interactions with the T-cell activation signalling pathway. Bioinformatics. (2015) 31:1075–83. doi: 10.1093/bioinformatics/btu787

PubMed Abstract | CrossRef Full Text | Google Scholar

137. Lindell DM, Berlin AA, Schaller MA, Lukacs NW. B cell antigen presentation promotes Th2 responses and immunopathology during chronic allergic lung disease. PLoS ONE. (2008) 3:e3129. doi: 10.1371/journal.pone.0003129

PubMed Abstract | CrossRef Full Text | Google Scholar

138. Zambrano S, De Toma I, Piffer A, Bianchi ME, Agresti A. NF-κB oscillations translate into functionally related patterns of gene expression. eLife. (2016) 5:E9100. doi: 10.7554/eLife.09100

PubMed Abstract | CrossRef Full Text | Google Scholar