The transitions leading from experimental data to the Boolean computational model are summarized in Fig. 1. Shown in Fig. 1 (Top) are four types of experimental data: spatial and temporal data on regulatory gene expression; results of trans perturbation experiments in which the activity of each regulatory gene is interrupted and the effects on other regulatory genes measured; and results of functional cis-regulatory analyses of specific interactions. In previous studies, the expression of every regulatory gene in the endomesoderm GRN model has been perturbed and the effects on expression of every other gene in the system measured (2–4). The experimental data give rise to two kinds of abstraction, both Boolean in form. Fig. 1 (Middle) shows that these are a matrix of spatial and temporal Boolean gene expression and a matrix of Boolean perturbation results. These matrices provide the nodes and interactions, which, by application of appropriate logic rules, are then integrated into the GRN model. The Boolean computational model we have constructed here represents a further level of abstraction, as shown in Fig. 1 (Bottom). Here we see that the predicted interactions at each gene in the GRN model were used to generate the logic driving the Boolean computational model. The principle used was as follows: the input regulatory information given in the GRN model at each node was used to formulate a logic equation that would express the activity of the gene as a function of its regulatory inputs. In the Boolean model, the regulatory state is thus the output of all logic equations at each point in time and in each spatial domain considered. Thus, were the GRN model essentially sufficient, the logic equations for all network nodes would ideally contain the necessary information to reproduce all observed regulatory states. As we discuss in the following sections, formalisms were developed for the computation of regulatory interactions between network nodes, for the relative positions of different embryonic domains, for the functions of transcription factors responding to signaling cascades, and for kinetic aspects of gene expression. The Boolean model computes expression of genes at every network node through time and space. The result permits a direct comparison between the observed and the computed Boolean gene expression profiles.

Open in a separate window

Fig. 1.

Diagram summarizing the flow of information and logic relations between experimental data, the GRN model deduced from it, and the Boolean computational model. Gene expression data, systemwide perturbation experiments, and cis-regulatory analyses (Top), combined with the basic facts of the embryonic process, are used to generate abstract time-space expression pattern and the interaction matrices (Middle). These results in turn underlie the design of the deduced genetic circuitry captured in the GRN model. Inputs to the Boolean computational model (Bottom) are the vector equations derived directly from each node of the GRN model; the relative geometry of the interacting embryonic spatial domains considered in the model; and the kinetics with which gene cascades operate in this embryo. The model incorporates essential aspects of the temporal and spatial biology of the embryo, and gives rise to a matrix of specific predictions of where, when, and for how long every gene is individually expressed over a 24-h period. This computed matrix of expression results can now be compared directly with the observed matrix of gene expressions, in normal or perturbed conditions.

Our object was to generate a stepwise computation in which, at every hour in real time, the transcriptional states—on or off—of all regulatory genes in the system are assessed, based on the computed outputs at the previous step. To this end, a computational application (GeNeTool) was developed. A flowchart describing the use of the logic equations in GeNeTool is shown in Fig. S1. Boolean computational models have been built to analyze other specific developmental processes, including expression patterns of segment polarity genes in the late cleavage Drosophila embryo (6, 7), and in mouse midhindbrain boundary formation (8). The model we describe here has many features that allow us to capture the causal impact of cis-regulatory function, signaling interactions, and cell lineage geometry, all of which contribute to the regulatory developmental biology of the sea urchin embryo. In brief, the model operates as follows: each node of the GRN model is computationally represented in the computational model by explicit statements, vector equations, capturing its cis-regulatory interactions. At each node, these statements are used to compute the nodal output at each step in time and in each spatial domain. These outputs are then fed to all the nodes in the system at the next step, and the computation is repeated. It is important to realize that this is a synthetic de novo computation and not a simulation in the sense that it simply restates in different language what is already included in the previous model. Thus, the starting GRN has no dynamics, nor does it encompass the geometrical features required for integration of signaling with gene expression as we describe later; nor is the output iteratively fit to the observed expression patterns as in a typical simulation. The goal here is not to simply reproduce observed expression patterns with a specific set of network interactions. The goal is instead to test the experimentally derived GRN topology and identify lacunae in its informational content. The outcome is comparable to an automaton that, when it has been fed the initial conditions, runs by itself until the end of the process considered.

The features of the developmental regulatory system served as direct guiding principles for the organization of our model, as follows.

The main function of the genomic regulatory system for development is specification of spatial patterns of gene expression. The first test for the usefulness of the Boolean computational model is therefore its ability to reproduce the correct regulatory state for each regulatory domain. The regulatory state was defined for each of the four endomesodermal regulatory domains as the sum of all regulatory genes expressed in that domain for at least one time interval between 18 and 30 h, except for the early mesoderm GRN, for which a 12- to 16-h time window was used because the later GRN model is not complete. The comparison between computed and observed spatial regulatory states is shown in Fig. 2. Transcriptional control is not yet understood for a few regulatory genes in the GRN model, and no vector equations could be formulated for these genes (Fig. 2, genes on white background). They were nonetheless incorporated into the Boolean model as they are known to provide inputs for other genes. For all J factors functioning downstream of signaling interactions, plus for genes predicted to be included in the GRN circuitry but not yet identified, there are no observed gene expression data available (Fig. 2, genes on black background). However, we have generated observed and computed gene expression data for 33 regulatory genes in the four domains, allowing 132 direct comparisons (Fig. 2, genes on gray background). Of these comparisons, all but two showed a perfect match between observed and computed spatial expression (Fig. 2). One exception was brn1/2/4, for which the computational model predicted expression in the early mesoderm, where expression of this gene has not been observed, in addition to the correctly computed expression in the veg2 endoderm. The reason is that all currently known regulatory inputs into brn1/2/4 are present in the early mesoderm, predicting that additional factor(s) must regulate brn1/2/4. The second deviation between computation and observation is caused by wnt8, which fails to clear from the veg2 endoderm domain in the computational model, indicating a missing repressor in the GRN model. These two exceptions aside, Fig. 2 proves that the regulatory information included in the model suffices to explain almost every known spatial regulatory gene expression pattern in the endomesoderm up to 30 h of development.

Open in a separate window

Fig. 2.

Comparison of computed and observed spatial expression patterns. Upper: All regulatory genes included in the model. Experimental expression data are summarized digitally online (http://sugp.caltech.edu/endomes/#BioTapestryViewer). Expression in this chart indicates presence of the computed or observed transcript at any time between 18 and 30 h, except for the V2 mesoderm domain, where the relevant interval is 12 to 16 h (see Computed Spatial Gene Expression above). As indicated in the key, the actual comparison between observed and computed results, i.e., for genes in which both types of data exist, pertains to the colored and gray squares. A discrepancy between computed and observed expression is indicated by a black bar. For genes or regulatory operators shown on black backgrounds, no observational data are available; for genes on white backgrounds, no upstream regulatory information is available, and the activation of these genes was not computed. Where such genes are expressed, or for maternal nuclear β-catenin, in which it is present, the time/space domain is shown in white; where they are observed not to be expressed is shown in gray. The disposition of the four spatial domains in the chart, in embryos viewed from the vegetal pole, is color-coded (Fig. S4 shows a general diagram of embryogenesis).

The comparison of computed and observed spatial expression of the many regulatory genes in Fig. 2 extends over a broad temporal window (18–30 h). However, when examined in detail, each gene displays its own individual temporal pattern of activation and sometimes repression. The temporal resolution of the observed expression patterns is based on whole mount in situ hybridization data, which exist at a time resolution of approximately 3 h for most genes included up to the 30 h developmental stage (data provided at http://sugp.caltech.edu/endomes/#BioTapestryViewer). In addition, a NanoString data set provides transcript accumulation measurements for all these genes at 1-h resolution (16).

The computed results are shown with respect to the observed expression patterns in Fig. 3, for every gene at every hour for every domain. Discrepancies between computed and observed expression in real time are considered significant if greater than the step time of 3 h (Fig. 3, solid black bars), whereas deviations less than 3 h (Fig. 3, open black bars) may not be real, as this is the limit of resolution of observed spatial data. We see, perhaps surprisingly, that the assumption of a uniform, canonical step time reproduces the observed dynamics within the experimental resolution for the vast majority of genes with no significant discrepancies from data (Fig. 3, colored and gray rectangles). In the skeletogenic domain, aside from a few 1-h discrepancies (beyond the 3-h observational variance), the model turned on foxb 3 h too early and failed to turn off foxn2/3 for 4 h when it should have; in the mesoderm domain, over the 6- to 18-h period, there was only one 1-h discrepancy; in the veg2 endoderm, brn1/2/4 was briefly turned on when it should not have been, and was later turned on a couple of hours too early, and tgif was turned on 3 h early; in veg1 endoderm, the only serious discrepancies were that otxβ was turned off 6 h too early and z13 was turned on for 4 h when it should not have been. For three genes, manual corrections were imposed to avoid promulgating downstream effects (Fig. 3, boxes with heavy lined black borders): delta and wnt8 cease to be expressed in veg2 endoderm after certain times for reasons not yet known, and these genes were turned off manually at the observed times; similarly, the known inputs into foxa (17) do not indicate why this gene continues to be transcribed after 24 h, and this activity was set manually for the 24- to 30-h interval. However, considering that there are 2,772 time/space gene expression domains computed in the model and shown in Fig. 3, this is a remarkably small frequency of significant temporal discrepancies, all of relatively brief duration, and there was only one spatial error.

Open in a separate window

Fig. 3.

Computed vs. observed spatial expression in real time. The genes are indicated at the top and the comparison symbolism is as in Fig. 2; i.e., plain filled colored or gray squares indicate exact congruence between computed and predicted expression. The four embryonic spatial domains (A–D) are separately portrayed, each against 1- to 30-h time scales. Open rectangles indicate discrepancies <3 h; as described in the text, as the observations have only 3-h resolution, these discrepancies may not be real. Discrepancies >3 h are indicated by black solid rectangles. The green striped bars indicate maternal mRNAs that are not part of the zygotic GRN model and therefore not computationally predicted. Some of these maternally supplied factors provide initial inputs to the computation. The three places where the whole column is striped (delta, foxa, and wnt8) indicate phases of expression or nonexpression known to occur, but not explained, and set manually in the computation.

The discrepancies are individually valuable items of information, for they pinpoint exactly what remains to be solved in terms of transcriptional causality. However, the main conclusion is that imposition of a uniform step time in the computation results in a generally accurate representation of the dynamics of the spatial expression patterns. The regulatory genes of the sea urchin embryo operate more or less similarly to one another, and, like a clock, the pace of which is set by the kinetics of the basic processes of transcription molecular biology. The application of a correct step time is important for the overall performance of the Boolean model. To test the implications of variations in the step time, we ran the model by using different step times. The catastrophic performance of the model when the step time was changed from 3 h to 4 h is shown in Fig. S7. Thus, we discovered that the embryo regulatory system operates largely like a clock by using a 3 h step time. This fact obviates detailed, gene-by-gene synthesis/turnover kinetic analysis. The few discrepancies most likely indicate missing regulatory inputs, rather than specific kinetic deviations from this general conclusion.

Sea urchin embryos have been experimentally exposed to various perturbations. These types of experiment offer an opportunity to challenge the completeness of the explanatory power of the Boolean model, again by comparing computational prediction to observed gene expression or phenotypic results, but now under specifically altered conditions that mimic in silico well studied experimental perturbations. Here we consider four such projects.

In the first of these, we mimicked extinction of delta expression by injection into the egg of delta morpholinos (4, 18, 19). The experimental result is loss of almost all mesodermal gene expression from the veg2 cells, which normally give rise to all nonskeletogenic mesoderm (i.e., the ring of cells shown in blue in Figs. 2 and ​and4).4). Instead, as symbolized in Fig. 4A (which represents the embryo at approximately 18 h), these cells express veg2 endodermal genes, which are now transcribed ectopically in a double rather than single ring of cells (4). When delta gene expression is manually turned off in the model, the computational result is precisely consistent with the experimental result (Fig. 4A; Fig. S8A shows complete detailed results). Only the two early mesoderm genes gcm and gatae have been examined in delta morpholino embryos, affording a direct comparison with data, and the expected loss of expression of these genes was obtained in silico. Similarly, although only hox11/13b, foxa, and blimp1 among veg2 endodermal genes have been studied experimentally in this regimen, the Boolean model again produced the expected result of predicted ectopic expression in the presumptive mesodermal ring of cells. The model also predicted that the same happens for the other veg2 endoderm genes in the model that have not yet been examined experimentally.

Open in a separate window

Fig. 4.

In silico perturbation experiments. Left: Experimental perturbation mimicked in the computation and diagram of its general result. The five panels containing gray and black columns (Right) indicate respectively the in silico perturbation and its computational effects in the four embryonic domains. The specific genes affected in the computation are shown below: “Observed & Computed” indicates genes with a comparable experimentally observed change in expression; “Additional Predictions” indicate genes that have not been experimentally studied in these perturbations. (A) Interference with skeletogenic micromere expression of delta and (B) global overexpression of Pmar1. The detailed results as computed for all genes in the model are shown in Fig. S8 A and B.

In a second perturbation, we mimicked experimentally induced global expression of the Pmar1 repressor by injection of pmar1 mRNA into the egg (2, 20, 21). The observed experimental result is dramatic: the whole embryo is converted into a ball of mesenchymal cells expressing skeletogenic genes, as indicated diagrammatically (Fig. 4B). This is also just what occurred in silico when pmar1 is manually activated in all domains from 5 h on: eight skeletogenic genes have been shown experimentally to be globally expressed in pmar1 embryos, and, in the computation, all eight were expressed in all domains (Fig. 4B). An additional prediction is that the same should be true of the remaining skeletogenic genes as well (Fig. S8B shows detailed results). Furthermore, the Boolean model predicted that expression of mesodermal genes, veg2 endoderm genes, and veg1 endoderm genes would be turned off. This result is certainly consistent with the powerful global mesenchymal phenotype (20).

In a third perturbation, we directly mimicked a recent molecular biology experiment (4) in which the quantitative and spatial effects on the whole GRN of blocking hox11/13b expression with hox11/13b morpholino was determined throughout the pregastrular period. This was done by manually extinguishing hox11/13b expression in veg2 and veg1 endoderm (Fig. S8C). With the exception of two discrepancies that were merely 1-h temporal deviations, every observed change in gene expression was generated computationally.

A fourth perturbation was more challenging. Here we attempted to use the Boolean model to determine whether the computed regulatory relationships suffice to predict the results of a famous blastomere transplantation experiment reported by Hoerstadius in 1935 (22) and, 58 y later, further explored by us by using a molecular marker (23). In this experiment, the four fourth cleavage skeletogenic micromeres are removed from a donor embryo and transplanted to the animal pole of an otherwise normal early cleavage embryo possessing its own set of vegetal micromeres. The result is that a complete second gut, plus a second set of nonskeletogenic mesoderm lineages (18), is induced to form from the cells at the animal pole of the embryo (Fig. 5A, diagram). Essentially, this means that the signals emitted from the transplanted skeletogenic micromeres elicit the encoded developmental GRNs, and this suffices to recreate all the regulatory and spatial relationships that normally generate the concentric regulatory states of the vegetal pole of the embryo, as shown diagrammatically in Fig. 2. Were we to mimic the micromere transplantation experiment in our model, would it regenerate in a computationally naive patch of cells the same set of concentric mesodermal and endodermal regulatory states as it generates (Figs. 2 and ​and3)3) at the normal vegetal end? We modeled these naive cells by endowing them with the global maternal factors, but without the localized nuclear β-catenin that normally occurs in veg2 and veg1 lineages (24, 25). The “transplanted” micromeres were allowed all their own normal specification functions, as these are autonomous with respect to the rest of the embryo. The surrounding naïve cells were considered as forming three domains (Fig. 5B), the immediately adjacent cells (ring 1), the cells surrounding these (ring 2), and the cells surrounding ring 2 (ring 3). We now asked whether the in silico regulatory states of ring 1 would equal those of the mesoderm, and those of rings 2 and 3 the regulatory states of veg2 and veg1 endoderm, respectively. Fig. 5C shows that this is exactly what occurred; the only discrepancies are two genes that are dependent on inputs from two genes for which vector equations are lacking (dac, myc). If we manually activate dac and myc expression, there are no discrepancies.

Open in a separate window

Fig. 5.

Computational recreation of the micromere transplantation experiment. (A) Diagrammatic summary of experiment in which a set of micromeres is transplanted to the animal pole of an eight-cell embryo; by the next cleavage, the endogenous set of micromeres has also appeared; and, ultimately, a complete second gut and mesoderm are induced to form. (B) Conceptual organization of animal pole cells surrounding the descendants of the transplanted micromeres, at ~200- to 400-cell stage, in similar geometric relation as at the vegetal pole. (C) Chart comparing computed regulatory states of ring 1 vs. V2 mesoderm, ring 2 vs. V2 endoderm, and ring 3 vs. V1 endoderm; symbols are as in Figs. 2 and ​and44 and analysis was carried out for the 22- to 30-h interval. The comparison between computed expression patterns in the normal vegetal region of the embryo and computed expression patterns surrounding the in silico transplanted micromeres is almost perfect.

In considering the import of these in silico perturbations, we note that, except for the hox11/13b case, none of the experimental perturbation results we sought to reproduce were used in any specific way to construct the linkages in the GRN models on which the vector equations of Fig. S3 were based. Thus, a model based on functional analysis in normal development correctly predicts independent experimental results in unnatural contexts: this provides a proof of principle that the genetic interaction system in the model operates essentially as in the embryo. Furthermore, in each case, the model extends our spotty understanding of the consequences of the perturbations, which is limited by the genes the experimentalists happen to have investigated, to all of the genes in the whole system, producing, for example, the additional predicted consequences indicated in Fig. 4.

Spatial gene expression is both the driver and product of developmental mechanisms. Here we show that observed embryonic regulatory gene expression in space and time can be recreated, systemwide, by reducing these mechanisms to the underlying gene regulatory interactions. We apply a Boolean interpretation to regulatory gene expression and regulatory gene interaction. This approach is consistent with the visibly apparent Boolean patterns of spatial gene expression, a fundamental property of embryonic development. Given cells at given times express given regulatory genes at detectable levels or do not significantly express them. Diverse mechanisms conspire to produce this basic property.

First, in sea urchin embryos, the transcriptional cascade kinetics themselves contribute in a crucial way. Our earlier analyses showed that the point at which a downstream gene is activated with respect to the upstream drivers long precedes the attainment of driver steady state (5), and this means that the kinetics of gene cascades in the embryo are relatively level-insensitive and insensitive to synthesis and decay rates. The next gene in a sequence goes on as soon as there is sufficient driver transcription factor in the system. In accord with this, measurements show that levels of regulatory gene transcripts in different individual batches of eggs in these polymorphic animals frequently vary by twofold or more (16): the developmental gene regulatory system cares only about on and off, and level control is sloppy. Of course, as our high-resolution time courses for all genes included in this model show (16), the expression levels change through time. However, our model is not aimed at explication of expression kinetics, but rather at explication of spatial domain of expression, and, for this, what is causal is the availability of activators and repressors. As we have seen, in this system, it takes approximately 3 h from the time a gene is detectably activated to the time its product can affect an immediate target gene. Thereafter, the levels of the gene product usually continue to increase within the same spatial domain without change in the set of affected target genes.

Second, the boundaries of Boolean gene expression territories are in embryonic development usually set by GRN subcircuits that use repression. A variety of subcircuits of different design are known that accomplish this function (10, 11). Although repression itself is usually modeled as a continuous equilibrium function, this is a severe oversimplification for animal cells because, following binding of the repressor, secondary changes in chromatin structure and composition follow that may produce irreversible repression states that can last for the life of the cell. Thus, per se, repression behaves in a Boolean manner in space and in time.

Third, even when there are graded upstream driver inputs, Boolean gene expression outputs occur downstream as a result of particular GRN subcircuit designs that function hysteretically (10). To take one example, in the vertebrate neural tube, a dynamically changing Sonic Hedgehog gradient input is interpreted by the encoded GRN to produce sharp stripes of regulatory gene expression (26). Thus, whereas the input is graded, the spatial regulatory output is Boolean. The majority of such mechanisms basically devolve from the input/output functions generated by the designs of a variety of particular GRN subcircuits, which can be said to perform “Booleanization” functions. In general, signaling interactions in development produce essentially Boolean regulatory outputs. Most known developmental signaling pathways function by use of a response system that, in cells receiving the signal, promote target gene expression, and, in cells not receiving it, prevent it, by repressing the same target genes. This is a feature captured directly in the structure of our model, as described earlier.

The aim of this computational exercise was to assess the sufficiency of a system-level causal explanation for development. This explanation, couched in the format of the GRN model, provided the crucial logic functions of the computational model. The success of the computational model in reproducing developmental gene expression, shows that our knowledge of the underlying GRNs is relatively complete. What we do not know is specifically circumscribed. Lacunae in the model fall into two classes. For the several “white genes” of Fig. 2, we lack sufficient experimental cis-regulatory and perturbation data on which to base vector equations that could explain their activation, but because their downstream effects are known, these genes are included subsequently in the model. Second, the model executes a few miscomputations, in which the control operations built in fail to explain why a given gene turns off at a certain time, or why it continues to be transcribed (Fig. 3). Identification of these gaps in the GRN model is one of the useful outcomes of the computation. In addition, as we have seen, the GeNeTool Boolean model lends itself readily to predictive exploration of cis- and trans-perturbations of the GRN.

The most important general outcome of this work is that an assumed mechanism, based exclusively on experimentally established control functions, suffices to explain specifically almost all observed regulatory gene expression in space and time. This model contributes a step-by-step mechanistic understanding of how the complicated GRN shown in Fig. S2 actually works. There remains little room for causal explanation of spatial regulatory gene expression by any other mechanism. The formalism presented here provides a general means of representing in silico the genomic regulatory system underlying animal body plan development and further applications may be anticipated, ranging from evolution to synthetic developmental biology.

A computational platform, named GeNeTool, was constructed to enable the model computation. The graphics and design functionalities of GeNeTool are to be described in detail elsewhere. However, Fig. S1 shows a flowchart that describes the overall process of the computation, and the program can be obtained from the authors on request. A brief summary of the rules of logic operation, which accommodate temporal and relative spatial relationships in the model, is presented in Fig. S9.