ENHANCERS IN DEVELOPMENT AND HUMAN DISEASE

Most of our understanding of the function of enhancers comes from developmental biology or studies of the genetics of human disease. Human traits typically have a hereditary component but demonstrate complex patterns of inheritance. Genome-wide association studies (GWASs) have been widely used to identify complex trait loci and have identified more than 25,000 single-nucleotide polymorphisms (SNPs) that are significantly associated with variation in more than 700 traits and diseases (74). Although GWASs can still explain only a small fraction of the phenotypic variance (50), the list of validated regulatory mutations responsible for heritable susceptibility to diseases is growing at a steady rate.

The significant role of regulatory variation in complex trait heritability is underscored by the finding that the vast majority of trait-associated SNPs are non-exonic (49) and occur within putative regulatory elements far more often than expected by chance (30, 52). This suggests that disruption of regulatory function is a common mechanism by which noncoding sequence variants contribute to human disease. When a regulatory variant is identified, it is often hypothesized that the variant disrupts a transcription factor (TF) binding site, creates a new binding site, or both. In a recently elucidated example, a GWAS showed that the common SNP rs339331 increases prostate cancer risk (68) (odds ratio = 1.22, p = 1.6 × 10⁻¹²). Huang et al. (32) dissected this locus and showed that the risk SNP allele TTTTATGAG is bound by HOXB13, while the protective allele TTTCATGAG is not bound by HOXB13. This particular TF, in combination with FOXA1 and AR, activates RFX6 and promotes cell migration and metastatic disease (32). Since approximately 50 SNPs are typically in tight linkage disequilibrium with each GWAS-associated variant, similar detailed experimentation will be required to identify causal variants within disease-associated loci, but only a small number of these loci have been studied in detail.

A long-standing problem encountered when attempting to generalize known binding site disruptions is that the biological consequences of variation in a specific binding site are strongly dependent on both cell type and the neighboring local sequence context, which defines the combinatorial TF interactions with cell-specific cofactors. Because most TF binding sites are short and degenerate, there are usually thousands of what appear to be very good binding sites in the genome, yet only a fraction of these are occupied in a given cell type (7). Mutation of a binding site will have a functional consequence only in an occupied site. Therefore, the combinatorial code that determines cell-specific TF occupancy will determine which variants can alter regulatory element activity. Several computational methods have shown promise in detecting and quantitatively assessing the impact of variants in enhancers by training on uniformly processed genome-wide epigenomic data sets generated by the Encyclopedia of DNA Elements (ENCODE) Consortium and National Institutes of Health Roadmap Epigenomics Mapping Consortium (20, 62, 77).

In the context of embryonic development, multicellular organisms require cells to make fate decisions by integrating extracellular cues ranging from biochemical to mechanical signals. We now have a sophisticated understanding of how extracellular inputs, especially signaling molecules such as SHH or TGF-β, are transduced intracellularly and ultimately activate a relatively small set of TFs that play pivotal roles in cell fate determination. For instance, Nodal/TGF-β signaling is required for specifying definitive endoderm (DE) differentiation in gastrulating mouse and zebrafish embryos (2, 14, 21, 63), a process that has been recapitulated using human embryonic stem cells (hESCs) and mouse embryonic stem cells through directed differentiation (17, 38). Nodal/TGF-β signaling activates the SMAD2–4 TFs, which then cooperate with key lineage TFs, including FOXH1, EOMES, MIXL1, and GATA6, to activate the DE transcriptional program (46, 48, 58, 69, 72).

Indeed, Li et al. (46) uncovered all of these core TF genes (FOXH1, EOMES, MIXL1, GATA6, SMAD2, and SMAD4), along with additional new regulators, in pooled, genome-scale CRISPR/Cas9 loss-of-function screens for genes that are required for DE differentiation from hESCs, as shown in Figure 1. This study used the ESC–DE differentiation system to interrogate DNA elements in the regulatory network that controls induction of hESCs into DE, induced in this system by TGF-β and Wnt signaling. DE differentiation was triggered when hESCs were treated with CHIR-99021 and activin A (Figure 1). A SOX17/GFP knock-in allele reported the DE fate (35), as assessed by flow cytometry. To detect regulators of this process, iCas9 SOX17/GFP cells were infected with the human Genome-Scale CRISPR Knock-Out (GeCKO) v2 guide RNA (gRNA) library. After selection for cells with viral integration and induction of Cas9 expression, DE differentiation was performed, and SOX17/GFP+ DE and SOX17/GFP− non-DE cells were isolated by fluorescence-activated cell sorting. The abundance of individual gRNAs in each population was determined by high-throughput sequencing: gRNAs that target positive or negative regulators of endoderm specification should be depleted or enriched, respectively, in SOX17/GFP+ compared with SOX17/GFP– cells. A Z-score was calculated for each gRNA based on the ratio of gRNA reads in the populations. The top 20 hits included almost all of the nonredundant, cell-autonomous required genes in the Nodal pathway (ACVR1B, SMAD2, and FOXH1) (63) as well as the established DE TF genes EOMES and MIXL1 (80). TFs required to maintain the ESC state can also be screened by sequencing the pool of gRNAs enriched in self-renewing conditions.

However, there is a major gap in our knowledge of how TFs control the gene regulatory networks that dictate cell fate decisions. We often know which TFs are required for the acquisition or maintenance of a cell state during development, but the exact cascade of molecular events that either drives the cell state transition or stabilizes the cell state is unclear. Genomic data such as chromatin immunoprecipitation sequencing (ChIP-seq) data can provide rich information regarding the chromatin association of a TF, and knockout studies can identify genes with altered expression levels when the TF is deleted. However, these experiments may not indicate direct transcriptional consequences. There are two challenges to establish the causality of the TF binding and gene expression changes relevant to cell fate determination. First, a TF (e.g., TF A) usually has multiple binding sites near a gene of interest (e.g., gene X). Thus, even if the deletion of TF A causes a change of gene X expression, the impact (if any) of individual TF A binding sites on the control of gene X expression is typically unknown. Second, differentiation involves a cascade of molecular events, so it is conceivable that TF A regulates TF B expression, which then directly regulates the expression of gene X, even though TF A may also bind to genomic regions near gene X. In fact, multiple TFs often bind to the same region, but they may or may not directly contribute to transcriptional regulation, and some of the TFs may have overlapping, additive, synergistic, or buffering effects. Globally, it is challenging to determine from the genomic occupancy pattern alone which TFs are required for the cell to make a specific cell fate decision (51).

Therefore, in order to establish a predictive gene regulatory network for cell fate control, it is necessary to build on our knowledge of the TFs required for fate specification. This will allow us to identify cognate functional enhancer regions and measure the local and global consequences of perturbing these enhancers. For this purpose, we have been focusing on enhancers that mediate the ESC–DE transition because of their importance to the development of endoderm-derived organs, including the pancreas and liver. We expect that the identification of functional TF binding sites within noncoding regulatory elements will establish edges (causal regulatory interactions between genes) in the gene regulatory networks. Measurement of the dynamics of these regulatory networks will form the basis for building quantitative models of enhancer function. Ultimately, these models should describe how combinatorial TF genomic occupancy at multiple enhancers controls lineage-specific gene expression in embryonic development, tissue homeostasis, regeneration, and aging.

CORE REGULATORY GENOMIC CIRCUITS

We believe it is useful to develop a conceptual framework to address the issues raised above. We reason that instead of characterizing or identifying enhancers based solely on their ability to drive gene expression, it would be more productive to devise targeted strategies to interrogate enhancers that play distinct roles in the gene regulatory network. In particular, some enhancers, by virtue of their ability to regulate the expression of fate-determining genes, may play central roles in development. Previous studies of developmental regulatory control have identified general principles that govern the organization and structure of gene regulatory networks, which are consistent across a wide range of multicellular model organisms (18). Particularly useful for our purpose is the separation of genes and enhancers into two classes: those whose primary role is to specify cell fate, which we will call core genes, and downstream genes, whose role lies downstream of the fate specification genes and which perform the necessary functions of the cell once its fate is set (Figure 2a). There are also therefore two classes of developmental enhancers based on their endogenous activity and position of the genes they regulate in the network: core enhancers, which have a global impact in terms of cell fate maintenance or transition, and peripheral enhancers, which regulate the expression of one or sometimes multiple adjacent genes in cis but have little or no global (developmental) impact. From the viewpoint of the gene regulatory networks that control developmental lineage decisions, the core enhancers are likely to regulate the expression of lineage-determining TFs that connect them into nonlinear networks that produce bifurcations between stable cell states. The peripheral enhancers are downstream targets that have an impact on specific gene expression levels (e.g., of differentiated cells) but do not feed back into the control circuit. This separation of core and peripheral regulators has been particularly useful in interpreting the heritability of complex traits (10, 13, 47).

The identification of separate core and peripheral genes and enhancers in the regulatory network has implications for how to model regulatory networks using computational genomics and machine learning and how to interrogate these classes of enhancers using distinct experimental methods. Computational machine learning and statistical methods rely on the existence of many examples that contain patterns that represent likely predictive causal mechanisms. In the case of enhancer prediction, the biological structure of gene regulatory circuits provides the redundant examples from which these patterns can be learned. Each peripheral gene is driven by a small set of enhancers (at least one per cell type in which the gene is expressed, and sometimes only one), each containing binding sites for the core TFs (18). The fact that core genes are greatly outnumbered by peripheral genes has two key consequences: First, the set of core TFs in each cell type (the TF vocabulary) is small enough that the TF vocabulary of a given cell type is of limited complexity and is learnable, and second, the large number of peripheral genes (say, 5,000–15,000) active in any cell type requires a large set of peripheral gene enhancers that can be used as examples to train the computational model. Thus, the success of machine learning methods in predicting peripheral enhancers by identification of core regulator binding sites within them (25, 42, 44) is consistent with the general principles derived from targeted studies of developmental regulatory networks (18).

Based on the study of the ESC–DE transition by Li et al. (46) and many previous works, we can summarize a schematic gene regulatory network model of the core regulators controlling the ESC–DE transition in Figure 2b. The features of this model are built upon and consistent with the following observations from both perturbative functional studies and computational sequence analysis of cell-specific distal enhancers in many cell types (as described in more detail below):

Examples of functional data across the ESC–DE transition supporting these global observations are shown in Figure 2c, using ChIP-seq and chromatin accessibility [assay for transposase-accessible chromatin using sequencing (ATAC-seq)] data at two core ESC regulator genes (OCT4 and SOX2) and two core DE regulator genes (SOX17 and GATA6).

The fact that DNA sequence–based modeling can accurately predict a held-out test set, coupled with the observation that the features required to make this classification map to a relatively small set of TFs, implies that the set of core regulatory TFs is small (comprising 5–20 TFs) and that the much larger set of enhancers containing these binding sites map to peripheral target genes that do not typically affect the activity of the core regulator TFs directly. Additionally, each target enhancer typically contains TF binding sites for multiple core regulator genes.

Classifying enhancers into core and peripheral groups enables them to be interrogated separately using different methods. For the core enhancers, one only needs to focus on perturbing putative enhancers that regulate a relatively small number of core TF genes. Because core enhancers regulate core lineage-determining TFs, their perturbation would have a global impact on cell fate decisions. It is therefore possible to use a downstream cell fate–specific reporter (e.g., OCT4/GFP for hESC state or SOX17/GFP for the DE state) and large-scale, pooled CRISPR perturbation screens to determine the impact of perturbation on the cell fate. The regulation of hESC self-renewal and hESC–DE transition may involve both overlapping and distinct cis-regulatory sequences. It should be feasible to identify cis-regulatory elements that regulate the reporter gene expression (OCT4 or SOX17). On the other hand, there are obvious risks for identifying enhancers that regulate an upstream lineage-determining TF, which in turn regulates the cell fate decision. This identification would require the enhancer to have not only a relatively large effect on the transcription of the target gene but also secondary, tertiary, or even relatively indirect effects that would ultimately affect the cell fate.

GENERATING ENHANCER SETS FROM EPIGENOMIC DATA

Computational machine learning and statistical methods rely on the existence of many data points or training examples from which to extract patterns that represent likely predictive causal mechanisms. In the case of enhancer prediction, the biological structure of these circuits provides the redundant training examples from which these patterns can be learned. To train a DNA sequence model of enhancer activity, an appropriate enhancer training set must first be generated.

Enhancer activity is associated with both increased chromatin accessibility and histone modifications to chromatin state that contribute to the establishment and maintenance of activity (4, 8, 16, 31, 60, 61). Many epigenomic functional assays interrogate this active state and generate peaks of activity that can be used to define putative enhancer sets to train enhancer models, including ATAC-seq (11), DNase I hypersensitive site sequencing (DNase-seq) (9, 15, 64, 65, 70), histone ChIP-seq for acetylation of histone H3 on lysine 27 (H3K27ac) or monomethylation of histone H3 on lysine 4 (H3K4me1), and TF ChIP-seq (73) when core TFs are known. Trimethylation of histone H3 on lysine 4 (H3K4me3) marks are typically specific to promoters. ATAC-seq and DNase-seq reflect chromatin accessibility, which does not necessarily indicate enhancer activity but does have the advantage of higher spatial resolution relative to histone marks, which tend to flank the core TF binding sites. In addition, one can use ATAC-seq and DNase-seq without knowing the complete set of relevant TFs. For a complete set of TF ChIP-seq experiments in a given cell type, the full complement of TFs must be known, and good antibodies must exist—a tall order.

Once a set of appropriate marks are chosen, there are two common approaches to defining the training set. Many methods train on a limited positive set of 10,000–20,000 cell-specific peaks (1, 3, 5, 22, 25, 29, 41–44, 53, 66, 78) spanning 100–1,000 base pairs centered on the peak, along with a negative set of equal size or larger. In our work, we have found that a set of 20,000 300-base-pair sequences is usually close to optimal. This cell-specific training set approach has the advantage of focusing on identifying the core TF binding sites in that cell type. Other methods (36, 79) bin the genome in regularly spaced fixed-length (1,000-base-pair) bins that are not necessarily centered on an activity peak but have a multiclass label reflecting the epigenomic state of that bin for the full set of training data sets (n = 919 for the DeepSEA framework). This regular training bin approach has the advantage of generating a large set of sequences required for deep neural networks (DNNs) to obtain strong class-label accuracy but may miss the subtleties in the differences in TF vocabulary between specific biologically relevant cell states. This may be a particular concern for the less well-covered cell types in ENCODE, and such models should be retrained for these cases. In particular, the RFX6 SNP example discussed above is missed by training on all ENCODE data sets at one time, even though the prostate cancer cell line LNCaP is included in the training set (6). When trained on focused ENCODE samples, sequence-based modeling can be used to refine the quality of the data sets (43).

DNA SEQUENCE–BASED MACHINE LEARNING ENHANCER MODELS

Support vector machines (SVMs) and DNNs are two of the main classes of machine learning methods that have been successful for enhancer prediction. These methods are trained to classify a set of positive and negative examples, and the result of training is a classifier score function that can make predictions for the class of sequences outside the training set. SVM and DNN methods differ in the way the classifier score function is specified and how the parameters of this function are determined from the training data. They also differ in how the DNA sequence nucleotides are converted into an input vector of mathematical feature scores for each sequence element to be classified. In gkm-SVM (25), for example, each sequence is converted to a normalized vector of integer gapped k-mer counts. The parameters of this k-mer vocabulary are specified before training. In our studies, we use the full list of gapped k-mers of length L with k informative positions and L − k free positions (gaps or wild cards) and typically use (L,k) = (10,6) or (11,7); the latter is slightly more accurate and approximately half as fast. In the DNN, the input feature is later usually converted into a 4 × L binary integer matrix, with each nucleotide represented by a permutation of (1,0,0,0).

The training sequence set determines the features detected and should be designed to most clearly reflect the specific biological processes one aims to model. For example, when building a sequence model to predict SNPs that affect chromatin accessibility [chromatin accessibility quantitative trait loci (caQTLs) or DNase I sensitivity quantitative trait loci (dsQTLs)] in a cell line or primary cells, it is important to include examples of all accessible regions that may be altered by genomic variants in the experiment. Thus, in order to predict dsQTLs in lymphoblasts (19, 42) or ATAC-QTLs in T cells (24), we trained gkm-SVM on a positive set of a large number (∼23,000) of peaks of length L = 300 base pairs centered on the peak signal versus a GC- and repeat-matched negative sequence set of the same size. For sets of this size, the gkm-SVM R package (27) is most convenient, but for larger training sets, LS-GKM (40) is recommended. Training the SVM yields a gkm-SVM score function specified by the set of support vector coefficients that optimally separate the sequence elements in the positive and negative training sets (25). The gapped k-mer weight distribution is constructed from the gapped k-mer counts in the support vectors, . We often map the gapped k-mer weights to full k-mer weights for ease of interpretability, which produces an equivalent scoring function after training (26). The tails of these weight distributions, shown in Figure 3a,b, encode the features required to distinguish the cell-specific enhancer activity in the positive and negative training sets. In this case, these weights encode the TF binding sites required to predict chromatin accessibility in lymphoblasts, and the long positive tail of the weight distribution from 1 to 6 in Figure 3b maps to binding sites for the 10 TFs shown in Figure 3c. In this case, there are three classes of features: CTCF, promoter-specific TF binding sites (NRF1, SP1, NFY, and ELK4), and lymphoblast distal enhancer-specific TF binding sites (IRF2, BATF, RUNX1, NF-κB, and PU.1). Lymphoblast dsQTL SNPs disrupt all of these TFs, and accurate prediction of lymphoblast dsQTLs requires resolution of all three classes of features. In DNNs, these important features are encoded in the first position weight matrix (PWM) layer of convolution filters (360 filters with a length of 8 base pairs are typically used).

An alternative training set design can be used to detect more specific regulatory signals. For instance, to detect the TFs controlling the ESC–DE transition, instead of training versus inaccessible genomic negative sequence and comparing the weight vectors, one can isolate the differentially active TFs by choosing as a positive set the most differentially accessible peaks. Figure 4a shows the ESC and DE ATAC-seq signal from the study by Li et al. (46) at the union of all peaks from each state. Training a gkm-SVM model using the 5,000 most differentially accessible peaks in DE as a positive set (blue in Figure 4a) and the 5,000 most differentially accessible peaks in an ESC as a negative set (red in Figure 4a) yields a classifier with an area under the receiver operating characteristics (AUROC) value of 0.92. The tails of this gkm-weight vector contain binding sites for the core TF regulators of the DE (TCF, EOMES, SMAD2/3, GATA, AP1, and FOXH1) and ESC states (OCT4, NANOG, SOX2, EBOX, and CTCF), whose PWMs and top weights are shown in Figure 4b,c.

Training on ATAC-seq data from human tissue–derived cells or differentiated stem cells often detects very similar regulatory programs, as shown for human islet and pancreatic progenitor cells in Figure 5. These islet-specific enhancers form islet-specific DNA looping interactions in a type 2 diabetes–associated locus, as shown in Figure 5a, as measured by promoter capture Hi-C (PCHi-C) (55). Although some computational models have been proposed to predict the gene targets of these enhancers (75), the predictive power of these methods is much lower than initially reported (12, 76), and additional higher-resolution enhancer–promoter interaction data are needed to develop improved models.

These sequence-based enhancer prediction models have been tested with luciferase and massively parallel reporter assays in a wide range of cell types (mouse liver, retina, neurons, and melanocytes as well as human T cells, lymphoblasts, and GM12878, K562, HepG2, and SK-N-SH cells) (6, 24, 34, 37, 39, 42, 53, 56, 59). Most recently, in a prediction assessment of massively parallel reporter assay data, Shigaki et al. (67) tested five enhancers and nine promoters using saturation mutagenesis in disease-relevant cell types and found that the best models combined sequence features derived from enhancer prediction models that were trained on different data sets from ENCODE and the Roadmap Epigenomics Mapping Consortium. Using the same approach as previous studies (25, 42), the authors trained gkm-SVM on only the cell type–relevant DNase I hypersensitive site and ATAC-seq data set, which produced an average overall correlation with expression output of 0.39, as shown in Figure 6a. Performance improved to a correlation of 0.58 when training on multiple data sets and combining the deltaSVM scores with random forest regression, as described in Reference 67. This method also improves correlations at individual enhancer and promoter loci, as shown in Figure 6b,c.

REGULATORY NETWORK MODELS OF CELL FATE TRANSITIONS

Li et al. (46) introduced a simple continuum model of a gene regulatory network for a bistable genetic switch that described some features of the ESC–DE transition using reaction rate equation models similar to those previously used to model cell state transitions (23, 33, 57). It is reasonable to question the form of these reaction rate equation models on theoretical grounds, since the continuum limit upon which they are based may not be completely valid for some TFs expressed at low levels. A more technically rigorous and computationally much more challenging approach would utilize stochastic Langevin (28) or chemical master equations (54). Nevertheless, properly modeled reaction rate equations yield much clearer interpretations that lead to more facile biological insights. Also, the agreement between the experimental perturbations and the initial modeling in the study by Li et al. (46), in addition to the precision and robustness of embryonic developmental cell state transitions and the stability of cell states critical to multicellular life, suggests that a more theoretically rigorous model would lead to qualitatively similar conclusions.

In the Li et al. (46) model, shown in Figure 7a,b, TF genes A and B activate their own transcription by binding to nearby enhancers that negatively regulate the other TF. We will use lowercase letters (a and b) to describe the genes and uppercase letters (A and B) for the protein products. Gene a is transcribed (Figure 7a) only when TF A is bound but TF B is not: When B is bound at the gene a locus, gene a is in a nonproductive transcriptional state, and thus the transcription rate of gene a is given by the bound concentration of A (in the absence of B) at its own gene, [aA]. We assume that the equilibrium occupancy of regulators A and B at genes a and b is established rapidly relative to rates of protein production, and therefore the equilibrium of A and B at their binding sites is given by Michaelis–Menten kinetics: [aA] = [a][A] / k_aA = [a₀][A] / (k_aA + [A] + (k_aA / k_aB)[B]), where k_aA is the dissociation constant for TF A at its gene a binding site, k_aB is the dissociation constant for TF B at its gene a binding site, and [a₀] is the total DNA concentration, which we will absorb into the transcription rate t_a. Similarly, [bB] is equal to [b₀][B] / (k_bB + [B] + (k_bB / k_bA)[A]). Both [A] and [B] are degraded at rate r, and the transcription of each is proportional to t_a[aA] and t_b[bB], yielding the model in Figure 7c.

Stochastic simulations of induced cell state transitions can be modeled by adding a time-dependent impulse of activin that increases the transcription of one TF (here, B is a DE-specific TF), and weakening auto-activation models the effect of JNK inhibition and allows transition to DE at lower activin concentrations, in agreement with experiment (46) (Figure 7d). The transcription rate for A is reduced by increased B, as shown in Figure 7e. For parameters t = 3.8, k = k_aA / k_aB = 2, this system equilibrates at either a high A–low B or high B–low A state, depending on the initial conditions, as shown in Figure 7f. The full stability of this model can be worked out simply when the parameters for A and B are symmetric, as in Figure 7g. In this case, there are four fixed points, as shown in Figure 7h, and both high x–low y and low x–high y states are stable if k = k_aA / k_aB > 1, in the usual case where t > 1. The fixed point at x = y is an unstable saddle for k > 1, so for k > 1 this system exhibits bistability and can transition from one state to another with a significant perturbation, as shown in Figure 7d and Reference 46. However, this stability analysis shows that this system is somewhat sensitive to parameter choices. Since k_aA and k_aB are the dissociation constants for A and B at gene a, k > 1 requires that the repressive TF B binds at gene a with stronger affinity than the activating TF A and vice versa at the gene b locus. While understandable in the context of this mathematical model, this seems to be a rather difficult and unnatural requirement to satisfy for all mammalian cellular circuits.

A more realistic model is shown in Figure 8, which now incorporates the observation from Figure 2 that there are multiple core TFs active in each cell state that bind cooperatively at activating and repressive regulatory DNA elements. Here, TFs A, B, and C bind cooperatively at their cognate enhancers with binding sites for TFs A, B, and C at genes a, b, and c and genes x, y, and z. However, at genes a, b, and c, the complex binding produces a transcriptionally productive DNA looping conformation, while at genes x, y, and z it produces a transcriptionally inactive conformation. Similarly, TFs X, Y, and Z bind cooperatively; activate genes x, y, and z; and produce a transcriptionally inactive DNA conformation when they bind at genes a, b, or c. Now the relevant kinetic parameters are the dissociation constants for the ABC and XYZ complexes at the relevant gene enhancers (e.g., k_aABC and k_aXYZ). The model equations for this situation are shown in Figure 8c, and when simulated, they can produce a transition from the high X, Y, and Z state to the high A, B, and C state, as shown in Figure 8d, with a sufficiently large perturbation. This system can also be studied with phase-plane analysis techniques if we assume x = A = B = C and y = X = Y = Z, essentially modeling the complex concentrations, as in Figure 8f, and the system is now bistable for all k, as shown in Figure 8e. Stability analysis shows that when t > 1.9, there are two fixed points at and , and since the Jacobian here is independent of k, these are stable nodes for all k, as shown in Figure 8g. Thus, this cooperative model produces multiple stable states over a much wider range of kinetic parameter choices. It is possible that the ubiquity of cooperative binding by multiple TFs at mammalian developmental enhancers may have evolved because the state-switching dynamics of these cooperative circuits are more robust to binding site strength parameters than less cooperative gene regulatory networks involving single TFs.

SUMMARY AND FUTURE ISSUES

DNA sequence–based enhancer prediction methods and perturbative and functional studies are complementary methods that can be used to investigate the genetic regulatory networks controlling cell states and transitions between them. We have shown that these computational and experimental studies have detected features that are broadly consistent with each other. Predictive sequence features typically map to a relatively small set of core TF regulators, and cell-specific enhancers contain binding sites for multiple core TF regulators. When machine learning–based models are trained on these sets of enhancers as discriminative classifiers, they can predict the impact of mutations in reporter assays with reasonable quantitative accuracy. When these sequence-based models are used to design continuum dynamical models of genetic networks, these models can describe transitions between cellular states. Further analysis of these network models can yield insights into the features of genetic networks and the types of nonlinear interactions that support the stability and transitions between differentiated cellular states.

Yet many aspects of the modeling can be dramatically improved. Currently, linear and nonlinear classifiers often yield comparable overall accuracy, although there is ample evidence that nonlinear interactions between TFs contribute to enhancer activity and promoter interactions. Learning nonlinear interactions requires more training data, and it is likely that improved models in reduced feature spaces will be needed to allow models to learn the parameters of the nonlinear interactions between TFs and regulatory elements with existing amounts of training data. Methods of training on integrated large data sets can provide targeted training sequence data to detect more subtle regulatory events and should yield improved models. Higher-resolution three-dimensional chromatin looping measurements will likely be required to more fully understand and model the regulatory element interactions mediating promoter activation and repression.

disclosure statement

The authors are not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.