Dynamics Of Gene Regulatory Networks Research Articles

Single-cell dissection of aggression in honeybee colonies.

Understanding how genotypic variation results in phenotypic variation is especially difficult for collective behaviour because group phenotypes arise from complex interactions among group members. A genome-wide association study identified hundreds of genes associated with colony-level variation in honeybee aggression, many of which also showed strong signals of positive selection, but the influence of these 'colony aggression genes' on brain function was unknown. Here we use single-cell (sc) transcriptomics and gene regulatory network (GRN) analyses to test the hypothesis that genetic variation for colony aggression influences individual differences in brain gene expression and/or gene regulation. We compared soldiers, which respond to territorial intrusion with stinging attacks, and foragers, which do not. Colony environment showed stronger influences on soldier-forager differences in brain gene regulation compared with brain gene expression. GRN plasticity was strongly associated with colony aggression, with larger differences in GRN dynamics detected between soldiers and foragers from more aggressive relative to less aggressive colonies. The regulatory dynamics of subnetworks composed of genes associated with colony aggression genes were more strongly correlated with each other across different cell types and brain regions relative to other genes, especially in brain regions involved with olfaction and vision and multimodal sensory integration, which are known to mediate bee aggression. These results show how group genetics can shape a collective phenotype by modulating individual brain gene regulatory network architecture.

Inference of cell type-specific gene regulatory networks on cell lineages from single cell omic datasets

Cell type-specific gene expression patterns are outputs of transcriptional gene regulatory networks (GRNs) that connect transcription factors and signaling proteins to target genes. Single-cell technologies such as single cell RNA-sequencing (scRNA-seq) and single cell Assay for Transposase-Accessible Chromatin using sequencing (scATAC-seq), can examine cell-type specific gene regulation at unprecedented detail. However, current approaches to infer cell type-specific GRNs are limited in their ability to integrate scRNA-seq and scATAC-seq measurements and to model network dynamics on a cell lineage. To address this challenge, we have developed single-cell Multi-Task Network Inference (scMTNI), a multi-task learning framework to infer the GRN for each cell type on a lineage from scRNA-seq and scATAC-seq data. Using simulated and real datasets, we show that scMTNI is a broadly applicable framework for linear and branching lineages that accurately infers GRN dynamics and identifies key regulators of fate transitions for diverse processes such as cellular reprogramming and differentiation.

Estimation of Ordinary Differential Equation Models for Gene Regulatory Networks Through Data Cloning.

Ordinary differential equations (ODEs) are widely used for elucidating dynamic processes in various fields. One of the applications of ODEs is to describe dynamics of gene regulatory networks (GRNs), which is a critical step in understanding disease mechanisms. However, estimation of ODE models for GRNs is challenging because of inflexibility of the model and noisy data with complex error structures such as heteroscedasticity, correlations between genes, and time dependency. In addition, either a likelihood or Bayesian approach is commonly used for estimation of ODE models, but both approaches have benefits and drawbacks in their own right. Data cloning is a maximum likelihood (ML) estimation method through the Bayesian framework. Since it works in the Bayesian framework, it is free from local optimum problems that are common drawbacks of ML methods. Also, its inference is invariant for the selection of prior distributions, which is a major issue in Bayesian methods. This study proposes an estimation method of ODE models for GRNs through data cloning. The proposed method is demonstrated through simulation and it is applied to real gene expression time-course data.

Sampled-data Control of Probabilistic Boolean Control Networks: A Deep Reinforcement Learning Approach

The rise of reinforcement learning (RL) has guided a new paradigm: unraveling the intervention strategies to control systems with unknown dynamics. Model-free RL provides an exhaustive framework to devise therapeutic methods to alter the regulatory dynamics of gene regulatory networks (GRNs). This paper presents an RL-based technique to control GRNs modeled as probabilistic Boolean control networks (PBCNs). In particular, a double deep-Q network (DDQN) approach is proposed to address the sampled-data control (SDC) problem of PBCNs, and optimal state feedback controllers are obtained, rendering the PBCNs stabilized at a given equilibrium point. Our approach is based on options, i.e., the temporal abstractions of control actions in the Markov decision processes (MDPs) framework. First, we define options and hierarchical options and give their properties. Then, we introduce multi-time models to compute the optimal policies leveraging the options framework. Furthermore, we present a DDQN algorithm: i) to concurrently design the feedback controller and the sampling period; ii) wherein the controller intelligently decides the sampled period to update the control actions under the SDC scheme. The presented method is model-free and offers scalability, thereby providing an efficient way to control large-scale PBCNs. Finally, we compare our control policy with state-of-the-art control techniques and validate the presented results.

M6A RNA Methylation-Based Epitranscriptomic Modifications in Plasticity-Related Genes via miR-124-C/EBPα-FTO-Transcriptional Axis in the Hippocampus of Learned Helplessness Rats.

Impaired synaptic plasticity has been linked to dynamic gene regulatory network changes. Recently, gene regulation has been introduced with the emerging concept of unique N6-methyladenosine (m6A)-based reversible transcript methylation. In this study, we tested whether m6A RNA methylation may potentially serve as a link between the stressful insults and altered expression of plasticity-related genes. Expression of plasticity genes Nr3c1, Creb1, Ntrk2; m6A-modifying enzymes Fto, methyltransferase like (Mettl)-3 and 14; DNA methylation enzymes Dnmt1, Dnmt3a; transcription factor C/ebp-α; and miRNA-124-3p were determined by quantitative polymerase chain reaction (qPCR) in the hippocampus of rats that showed susceptibility to develop stress-induced depression (learned helplessness). M6A methylation of plasticity-related genes was determined following m6A mRNA immunoprecipitation. Chromatin immunoprecipitation was used to examine the endogenous binding of C/EBP-α to the Fto promoter. MiR-124-mediated post-transcriptional inhibition of Fto via C/EBPα was determined using an in vitro model. Hippocampus of learned helplessness rats showed downregulation of Nr3c1, Creb1, and Ntrk2 along with enrichment in their m6A methylation. A downregulation in demethylating enzyme Fto and upregulation in methylating enzyme Mettl3 were also noted. The Fto promoter was hypomethylated due to the lower expression of Dnmt1 and Dnmt3a. At the same time, there was a lower occupancy of transcription factor C/EBPα on the Fto promoter. Conversely, C/ebp-α transcript was downregulated via induced miR-124-3p expression. Our study mechanistically linked defective C/EBP-α-FTO-axis, epigenetically influenced by induced expression of miR-124-3p, in modifying m6A enrichment in plasticity-related genes. This could potentially be linked with abnormal neuronal plasticity in depression.

Deciphering the Transcriptional Landscape of Human Pluripotent Stem Cell-Derived GnRH Neurons: The Role of Wnt Signaling in Patterning the Neural Fate.

Hypothalamic gonadotropin-releasing hormone (GnRH) neurons lay the foundation for human development and reproduction; however, the critical cell populations and the entangled mechanisms underlying the development of human GnRH neurons remain poorly understood. Here, by using our established human pluripotent stem cell-derived GnRH neuron model, we decoded the cellular heterogeneity and differentiation trajectories at the single-cell level. We found that a glutamatergic neuron population, which generated together with GnRH neurons, showed similar transcriptomic properties with olfactory sensory neuron and provided the migratory path for GnRH neurons. Through trajectory analysis, we identified a specific gene module activated along the GnRH neuron differentiation lineage, and we examined one of the transcription factors, DLX5, expression in human fetal GnRH neurons. Furthermore, we found that Wnt inhibition could increase DLX5 expression and improve the GnRH neuron differentiation efficiency through promoting neurogenesis and switching the differentiation fates of neural progenitors into glutamatergic neurons/GnRH neurons. Our research comprehensively reveals the dynamic cell population transition and gene regulatory network during GnRH neuron differentiation.

The flowering transition pathways converge into a complex gene regulatory network that underlies the phase changes of the shoot apical meristem in Arabidopsis thaliana.

Post-embryonic plant development is characterized by a period of vegetative growth during which a combination of intrinsic and extrinsic signals triggers the transition to the reproductive phase. To understand how different flowering inducing and repressing signals are associated with phase transitions of the Shoot Apical Meristem (SAM), we incorporated available data into a dynamic gene regulatory network model for Arabidopsis thaliana. This Flowering Transition Gene Regulatory Network (FT-GRN) formally constitutes a dynamic system-level mechanism based on more than three decades of experimental data on flowering. We provide novel experimental data on the regulatory interactions of one of its twenty-three components: a MADS-box transcription factor XAANTAL2 (XAL2). These data complement the information regarding flowering transition under short days and provides an example of the type of questions that can be addressed by the FT-GRN. The resulting FT-GRN is highly connected and integrates developmental, hormonal, and environmental signals that affect developmental transitions at the SAM. The FT-GRN is a dynamic multi-stable Boolean system, with 223 possible initial states, yet it converges into only 32 attractors. The latter are coherent with the expression profiles of the FT-GRN components that have been experimentally described for the developmental stages of the SAM. Furthermore, the attractors are also highly robust to initial states and to simulated perturbations of the interaction functions. The model recovered the meristem phenotypes of previously described single mutants. We also analyzed the attractors landscape that emerges from the postulated FT-GRN, uncovering which set of signals or components are critical for reproductive competence and the time-order transitions observed in the SAM. Finally, in the context of such GRN, the role of XAL2 under short-day conditions could be understood. Therefore, this model constitutes a robust biological module and the first multi-stable, dynamical systems biology mechanism that integrates the genetic flowering pathways to explain SAM phase transitions.

Dynamic regulatory module networks for inference of cell type-specific transcriptional networks.

Changes in transcriptional regulatory networks can significantly alter cell fate. To gain insight into transcriptional dynamics, several studies have profiled bulk multi-omic data sets with parallel transcriptomic and epigenomic measurements at different stages of a developmental process. However, integrating these data to infer cell type–specific regulatory networks is a major challenge. We present dynamic regulatory module networks (DRMNs), a novel approach to infer cell type–specific cis-regulatory networks and their dynamics. DRMN integrates expression, chromatin state, and accessibility to predict cis-regulators of context-specific expression, where context can be cell type, developmental stage, or time point, and uses multitask learning to capture network dynamics across linearly and hierarchically related contexts. We applied DRMNs to study regulatory network dynamics in three developmental processes, each showing different temporal relationships and measuring a different combination of regulatory genomic data sets: cellular reprogramming, liver dedifferentiation, and forward differentiation. DRMN identified known and novel regulators driving cell type–specific expression patterns, showing its broad applicability to examine dynamics of gene regulatory networks from linearly and hierarchically related multi-omic data sets.

Kinetic networks identify key regulatory nodes and transcription factor functions in early adipogenesis

Sequence‐specific transcription factors (TFs) bind DNA, modulate chromatin structure, regulate gene expression, and orchestrate transcription cascades. Activation and repression of TFs drive tightly controlled regulatory programs that lead to cellular processes such as differentiation. We measured chromatin accessibility and nascent transcription at seven time points over the first four hours of induced adipogenesis of preadipocytes to construct dynamic gene regulatory networks. Regulatory networks describe successive waves of TF binding and dissociation followed by direct regulation of proximal genes. We identified 14 families of TFs that coordinate with and antagonize each other to regulate early adipogenesis. We developed a compartment model to quantify individual TF contributions to RNA polymerase initiation and pause release rates. Network analysis showed that the glucocorticoid receptor and AP1 drive immediate gene activation, including induction of Twist2. Twist2 is a highly interconnected node within the network and its expression leads to repression of target genes. Although Twist2's role in adipogenesis has not been previously appreciated, both Twist2 knockout mice and Setleis syndrome (Twist2‐/‐) patients lack subcutaneous and brown adipose tissue. We found that kinetic networks integrating chromatin structure and nascent transcription dynamics identify key genes, TF functions, and coordinate interactions within regulatory cascades.

Reconstruction of dynamic regulatory networks reveals signaling-induced topology changes associated with germ layer specification.

SummaryElucidating regulatory relationships between transcription factors (TFs) and target genes is fundamental to understanding how cells control their identity and behavior. Unfortunately, existing computational gene regulatory network (GRN) reconstruction methods are imprecise, computationally burdensome, and fail to reveal dynamic regulatory topologies. Here, we present Epoch, a reconstruction tool that uses single-cell transcriptomics to accurately infer dynamic networks. We apply Epoch to identify the dynamic networks underpinning directed differentiation of mouse embryonic stem cells (ESCs) guided by multiple signaling pathways, and we demonstrate that modulating these pathways drives topological changes that bias cell fate potential. We also find that Peg3 rewires the pluripotency network to favor mesoderm specification. By integrating signaling pathways with GRNs, we trace how Wnt activation and PI3K suppression govern mesoderm and endoderm specification, respectively. Finally, we identify regulatory circuits of patterning and axis formation that distinguish in vitro and in vivo mesoderm specification.

Data-driven modeling predicts gene regulatory network dynamics during the differentiation of multipotential hematopoietic progenitors.

Cellular differentiation during hematopoiesis is guided by gene regulatory networks (GRNs) comprising transcription factors (TFs) and the effectors of cytokine signaling. Based largely on analyses conducted at steady state, these GRNs are thought to be organized as a hierarchy of bistable switches, with antagonism between Gata1 and PU.1 driving red- and white-blood cell differentiation. Here, we utilize transient gene expression patterns to infer the genetic architecture-the type and strength of regulatory interconnections-and dynamics of a twelve-gene GRN including key TFs and cytokine receptors. We trained gene circuits, dynamical models that learn genetic architecture, on high temporal-resolution gene-expression data from the differentiation of an inducible cell line into erythrocytes and neutrophils. The model is able to predict the consequences of gene knockout, knockdown, and overexpression experiments and the inferred interconnections are largely consistent with prior empirical evidence. The inferred genetic architecture is densely interconnected rather than hierarchical, featuring extensive cross-antagonism between genes from alternative lineages and positive feedback from cytokine receptors. The analysis of the dynamics of gene regulation in the model reveals that PU.1 is one of the last genes to be upregulated in neutrophil conditions and that the upregulation of PU.1 and other neutrophil genes is driven by Cebpa and Gfi1 instead. This model inference is confirmed in an independent single-cell RNA-Seq dataset from mouse bone marrow in which Cebpa and Gfi1 expression precedes the neutrophil-specific upregulation of PU.1 during differentiation. These results demonstrate that full PU.1 upregulation during neutrophil development involves regulatory influences extrinsic to the Gata1-PU.1 bistable switch. Furthermore, although there is extensive cross-antagonism between erythroid and neutrophil genes, it does not have a hierarchical structure. More generally, we show that the combination of high-resolution time series data and data-driven dynamical modeling can uncover the dynamics and causality of developmental events that might otherwise be obscured.

The identifiability of gene regulatory networks: the role of observation data.

Identifying gene regulatory networks (GRN) from observation data is significant to understand biological systems. Conventional studies focus on improving the performance of identification algorithms. However, besides algorithm performance, the GRN identification is strongly depended on the observation data. In this work, for three GRN S-system models, three observation data collection schemes are used to perform the identifiability test procedure. A modified genetic algorithm-particle swarm optimization algorithm is proposed to implement this task, including the multi-level mutation operation and velocity limitation strategy. The results show that, in scheme 1 (starting from a special initial condition), the GRN systems are of identifiability using the sufficient transient observation data. In scheme 2, the observation data are short of sufficient system dynamic. The GRN systems are not of identifiability even though the state trajectories can be reproduced. As a special case of scheme 2, i.e., the steady-state observation data, the equilibrium point analysis is given to explain why it is infeasible for GRN identification. In schemes 1 and 2, the observation data are obtained from zero-input GRN systems, which will evolve to the steady state at last. The sufficient transient observation data in scheme 1 can be obtained by changing the experimental conditions. Additionally, the valid observation data can be also obtained by means of adding impulse excitation signal into GRN systems (scheme 3). Consequently, the GRN systems are identifiable using scheme 3. Owing to its universality and simplicity, these results provide a guide for biologists to collect valid observation data for identifying GRNs and to further understand GRN dynamics.

Transcriptomic, proteomic, and phosphoproteomic analyses reveal dynamic signaling networks influencing long-grain rice development

The LGS1 (Large grain size 1) gene, also known as GS2/GL2/OsGRF4, is involved in regulating grain size and quality in rice, but the mechanism governing grain size has not been elucidated. We performed transcriptomic, proteomic, and phosphoproteomic analyses of young rice panicles in Samba (a wild-type cultivar with extra-small grain) and NIL-LGS1 (a nearly isogenic line of LGS1 with large grain in the Samba genetic background) at three developmental stages (4–6) to identify internal dynamic functional networks determining grain size that are mediated by LGS1. Differentially expressed proteins formed seven highly functionally correlated clusters. The concordant regulation of multiple functional clusters may be key features of the development of grain length in rice. In stage 5, 16 and 24 phosphorylated proteins were significantly up-regulated and down-regulated, and dynamic phosphorylation events may play accessory roles in determining rice grain size by participating in protein–protein interaction networks. Transcriptomic analysis in stage 5 showed that differentially expressed alternative splicing events and dynamic gene regulatory networks based on 39 transcription factors and their highly correlated target genes might contribute to rice grain development. Integrative multilevel omics analysis suggested that the regulatory network at the transcriptional and posttranscriptional levels could be directly manifested at the translational level, and this analysis also suggested a regulatory mechanism, regulation of protein translation levels, in the biological process that extends from transcript to protein to the development of grain. Functional analysis suggested that biological processes including MAPK signaling, calcium signaling, cell proliferation, cell wall, energy metabolism, hormone pathway, and ubiquitin-proteasome pathway might be involved in LGS1-mediated regulation of grain length. Thus, LGS1-mediated regulation of grain size is affected by dynamic transcriptional, posttranscriptional, translational and posttranslational changes.

Macroscopic dynamics of gene regulatory networks revealed by individual entropy

Complex systems are usually high-dimensional with intricate interactions among internal components, and may display complicated dynamics under different conditions. While it is difficult to measure detailed dynamics of each component, proper macroscopic description of a complex system is crucial for quantitative studies. In biological systems, each cell is a complex system containing a huge number of molecular components that are interconnected with each other through intricate molecular interaction networks. Here, we consider gene regulatory networks in a cell, and introduce individual entropy as a macroscopic variable to quantify the transcriptional dynamics in response to changes in random perturbations and/or network structures. The proposed individual entropy measures the information entropy of a system at each instant with respect to a basal reference state, and may provide temporal dynamics to characterize switches of system states. Individual entropy provides a method to quantify the stationary macroscopic dynamics of a gene set that is dependent on the gene regulation connections, and can be served as an indicator for the evolution of network structure variation. Moreover, the individual entropy with reference to a preceding state enables us to characterize different dynamic patterns generated from varying network structures. Our results show that the proposed individual entropy can be a valuable macroscopic variable of complex systems in characterizing the transition processes from order to disorder dynamics, and to identify the critical events during the transition process.

Gene Regulatory Network Dynamical Logical Models for Plant Development.

Mathematical and computational approaches that integrate and model the concerted action of multiple genetic and nongenetic components holding highly nonlinear interactions are fundamental for the study of developmental processes. Among these, gene regulatory network (GRN) dynamical models are very useful to understand how diverse types of regulatory constraints restrict the multigene expression patterns that characterize different cell fates. In this chapter we present a hands-on approach to model GRN dynamics, taking as a working example a well-curated and experimentally grounded GRN developmental module proposed by our group: the flower organ specification gene regulatory network (FOS-GRN). We demonstrate how to build and analyze a GRN model according to the following steps: (1) integration of molecular genetic data and formulation of logical rules specifying the dynamic behavior of each gene; (2) determination of steady states (attractors) corresponding to each cell type; (3) validation of the GRN model; and (4) extension of the deterministic model with the inclusion of stochasticity in order to model cell-state transitions dependent on noise due to fluctuations of the involved gen products. The methodologies explained here in detail can be applied to any other developmental module.

Optimal control of probabilistic Boolean control networks: A scalable infinite horizon approach

AbstractOne of the major issues in systems biology is developing control theory for gene regulatory networks (GRNs). Particularly, an important objective is to develop therapeutic intervention strategies to alter regulatory dynamics of GRNs to avoid undesired or diseased states. Several optimal control strategies have been developed to find therapeutic intervention strategies for small (or medium) sized GRNs modeled as probabilistic Boolean control networks (PBCNs). However, due to the humongous nature of GRNs, we require an optimal control strategy that scales to large systems without posing any constraints on the network dynamics. In this article, we formulate the infinite horizon discounted cost optimal control problem by leveraging the Markov decision process (MDP) based PBCN structure to model GRNs. Further, we design a stationary optimal control strategy, thereby avoiding the states associated with a disease. An augmented state space is defined to facilitate the control design method. By exploiting Kullback–Leibler (KL) divergence and the log transformation of the cost function referred to as the desirability function, the optimal control is devised using the path integral (PI) approach. We propose a sampling‐based technique for the approximation of the PI and hence the optimal control of PBCNs. The sampling‐based strategy is amicable to parallel implementation, thereby addressing the optimal control problem of large‐scale PBCNs. Finally, the viability of the results in the article is shown by some illustrative examples.

Recombination facilitates genetic assimilation of new traits in gene regulatory networks.

A new phenotypic variant may appear first in organisms through plasticity, that is, as a response to an environmental signal or other nongenetic perturbation. If such trait is beneficial, selection may increase the frequency of alleles that enable and facilitate its development. Thus, genes may take control of such traits, decreasing dependence on nongenetic disturbances, in a process called genetic assimilation. Despite an increasing amount of empirical studies supporting genetic assimilation, its significance is still controversial. Whether genetic assimilation is widespread depends, to a great extent, on how easily mutation and recombination reduce the trait's dependence on nongenetic perturbations. Previous research suggests that this is the case for mutations. Here we use simulations of gene regulatory network dynamics to address this issue with respect to recombination. We find that recombinant offspring of parents that produce a new phenotype through plasticity are more likely to produce the same phenotype without requiring any perturbation. They are also prone to preserve the ability to produce that phenotype after genetic and nongenetic perturbations. Our work also suggests that ancestral plasticity can play an important role for setting the course that evolution takes. In sum, our results indicate that the manner in which phenotypic variation maps unto genetic variation facilitates evolution through genetic assimilation in gene regulatory networks. Thus, we contend that the importance of this evolutionary mechanism should not be easily neglected.

Scalable control of developmental timetables by epigenetic switching networks.

During development, progenitor cells follow timetables for differentiation that span many cell generations. These developmental timetables are robustly encoded by the embryo, yet scalably adjustable by evolution, facilitating variation in organism size and form. Epigenetic switches, involving rate-limiting activation steps at regulatory gene loci, control gene activation timing in diverse contexts, and could profoundly impact the dynamics of gene regulatory networks controlling developmental lineage specification. Here, we develop a mathematical framework to model regulatory networks with genes controlled by epigenetic switches. Using this framework, we show that such epigenetic switching networks uphold developmental timetables that robustly span many cell generations, and enable the generation of differentiated cells in precisely defined numbers and fractions. Changes to epigenetic switching networks can readily alter the timing of developmental events within a timetable, or alter the overall speed at which timetables unfold, enabling scalable control over differentiated population sizes. With their robust, yet flexibly adjustable nature, epigenetic switching networks could represent central targets on which evolution acts to manufacture diversity in organism size and form.

Multiple morphogens and rapid elongation promote segmental patterning during development.

The vertebrate hindbrain is segmented into rhombomeres (r) initially defined by distinct domains of gene expression. Previous studies have shown that noise-induced gene regulation and cell sorting are critical for the sharpening of rhombomere boundaries, which start out rough in the forming neural plate (NP) and sharpen over time. However, the mechanisms controlling simultaneous formation of multiple rhombomeres and accuracy in their sizes are unclear. We have developed a stochastic multiscale cell-based model that explicitly incorporates dynamic morphogenetic changes (i.e. convergent-extension of the NP), multiple morphogens, and gene regulatory networks to investigate the formation of rhombomeres and their corresponding boundaries in the zebrafish hindbrain. During pattern initiation, the short-range signal, fibroblast growth factor (FGF), works together with the longer-range morphogen, retinoic acid (RA), to specify all of these boundaries and maintain accurately sized segments with sharp boundaries. At later stages of patterning, we show a nonlinear change in the shape of rhombomeres with rapid left-right narrowing of the NP followed by slower dynamics. Rapid initial convergence improves boundary sharpness and segment size by regulating cell sorting and cell fate both independently and coordinately. Overall, multiple morphogens and tissue dynamics synergize to regulate the sizes and boundaries of multiple segments during development.

Temporal dynamics of a CSF1R signaling gene regulatory network involved in epilepsy.

Colony Stimulating Factor 1 Receptor (CSF1R) is a potential target for anti-epileptic drugs. However, inhibition of CSF1R is not well tolerated by patients, thereby prompting the need for alternative targets. To develop a framework for identification of such alternatives, we here develop a mathematical model of a pro-inflammatory gene regulatory network (GRN) involved in epilepsy and centered around CSF1R. This GRN comprises validated transcriptional and post-transcriptional regulations involving STAT1, STAT3, NFκB, IL6R, CSF3R, IRF8, PU1, C/EBPα, TNFR1, CSF1 and CSF1R. The model was calibrated on mRNA levels of all GRN components in lipopolysaccharide (LPS)-treated mouse microglial BV-2 cells, and allowed to predict that STAT1 and STAT3 have the strongest impact on the expression of the other GRN components. Microglial BV-2 cells were selected because, the modules from which the GRN was deduced are enriched for microglial marker genes. The function of STAT1 and STAT3 in the GRN was experimentally validated in BV-2 cells. Further, in silico analysis of the GRN dynamics predicted that a pro-inflammatory stimulus can induce irreversible bistability whereby the expression level of GRN components occurs as two distinct states. The irreversibility of the switch may enforce the need for chronic inhibition of the CSF1R GRN in order to achieve therapeutic benefit. The cell-to-cell heterogeneity driven by the bistability may cause variable therapeutic response. In conclusion, our modeling approach uncovered a GRN controlling CSF1R that is predominantly regulated by STAT1 and STAT3. Irreversible inflammation-induced bistability and cell-to-cell heterogeneity of the GRN provide a theoretical foundation to the need for chronic GRN control and the limited potential for disease modification via inhibition of CSF1R.