Body plan evolvability: The role of variability in gene regulatory networks.

Recent computational modeling of early fruit fly (Drosophila) development has characterized the degree to which gene regulation networks can be robust to natural variability. In the first few hours of development, broad spatial gradients of maternally derived transcription factors activate embryonic gap genes. These gap patterns determine the subsequent segmented insect body plan through pair-rule gene expression. Gap genes are expressed with greater spatial precision than the maternal patterns. Computational modeling of the gap-gap regulatory interactions provides a mechanistic understanding for this robustness to maternal variability in wild-type (WT) patterning. A long-standing question in evolutionary biology has been how a system which is robust, such as the developmental program creating any particular species' body plan, is also evolvable, i.e. how can a system evolve or speciate, if the WT form is strongly buffered and protected? In the present work, we use the WT model to explore the breakdown of such Waddington-type 'canalization'. What levels of variability will push the system out of the WT form; are there particular pathways in the gene regulatory mechanism which are more susceptible to losing the WT form; and when robustness is lost, what types of forms are most likely to occur (i.e. what forms lie near the WT)? Manipulating maternal effects in several different pathways, we find a common gap 'peak-to-step' pattern transition in the loss of WT. We discuss these results in terms of the evolvability of insect segmentation, and in terms of experimental perturbations and mutations which could test the model predictions. We conclude by discussing the prospects for using continuum models of pattern dynamics to investigate a wider range of evo-devo problems.

Evolution of the RNA world: From signals to codes.

The accumulated material in evolutionary biology, greatly enhanced by the achievements of modern synthetic biology, allows us to envision certain key hypothetical stages of prebiotic (chemical) evolution. This is often understood as the further evolution in the RNA World towards the RNA-protein World. It is a path towards the emergence of translation and the genetic code (I), signaling pathways with signaling molecules (II), and the appearance of RNA-based components of future gene regulatory networks (III). We believe that these evolutionary paths can be constructively viewed from the perspective of the concept of biological codes (Barbieri, 2003). Crucial evolutionary events in these directions would involve the emergence of RNA-based adaptors. Such adaptors connect two families of functionally and chemically distinct molecules into one functional entity. The emergence of primitive translation processes is undoubtedly the major milestone in the evolutionary path towards modern life. The key aspect here is the appearance of adaptors between amino acids and their cognate triplet codons. The initial steps are believed to involve the emergence of proto-transfer RNAs capable of self-aminoacylation. The second significant evolutionary breakthrough is the development of biochemical regulatory networks based on signaling molecules of the RNA World (ribonucleotides and their derivatives), as well as receptors and effectors (riboswitches) for these messengers. Some authors refer to this as the "lost language of the RNA World." The third evolutionary step is the emergence of signal sequences for ribozymes on the molecules of their RNA targets. This level of regulation in the RNA World is comparable to the gene regulatory networks of modern organisms. We believe that the signal sequences on target molecules have been rediscovered and developed by evolution into the gene regulatory networks of modern cells. In conclusion, the immense diversity of modern biological codes, in some of its key characteristics, can be traced back to the achievements of prebiotic evolution.

Problem of Domain/Building Block Preservation in the Evolution of Biological Macromolecules and Evolutionary Computation.

Structurally and functionally isolated domains in biological macromolecular evolution, both natural and artificial, are largely similar to "schemata", building blocks (BBs), in evolutionary computation (EC). The problem of preserving in subsequent evolutionary searches the already found domains / BBs is well known and quite relevant in biology as well as in EC. Both biology and EC are seeing parallel and independent development of several approaches to identifying and preserving previously identified domains / BBs. First, we notice the similarity of DNA shuffling methods in synthetic biology and multi-parent recombination algorithms in EC. Furthermore, approaches to computer identification of domains in proteins that are being developed in biology can be aligned with BB identification methods in EC. Finally, approaches to chimeric protein libraries optimization in biology can be compared to evolutionary search methods based on probabilistic models in EC. We propose to validate the prospects of mutual exchange of ideas and transfer of algorithms and approaches between evolutionary systems biology and EC in these three principal directions. A crucial aim of this transfer is the design of new advanced experimental techniques capable of solving more complex problems of in vitro evolution.

Heuristic algorithms in evolutionary computation and modular organization of biological macromolecules: Applications to in vitro evolution.

Evolutionary computing (EC) is an area of computer sciences and applied mathematics covering heuristic optimization algorithms inspired by evolution in Nature. EC extensively study all the variety of methods which were originally based on the principles of selectionism. As a result, many new algorithms and approaches, significantly more efficient than classical selectionist schemes, were found. This is especially true for some families of special problems. There are strong arguments to believe that EC approaches are quite suitable for modeling and numerical analysis of those methods of synthetic biology and biotechnology that are known as in vitro evolution. Therefore, it is natural to expect that the new algorithms and approaches developed in EC can be effectively applied in experiments on the directed evolution of biological macromolecules. According to the John Holland's Schema theorem, the effective evolutionary search in genetic algorithms (GA) is provided by identifying short schemata of high fitness which in the further search recombine into the larger building blocks (BBs) with higher and higher fitness. The multimodularity of functional biological macromolecules and the preservation of already found modules in the evolutionary search have a clear analogy with the BBs in EC. It seems reasonable to try to transfer and introduce the methods of EC, preserving BBs and essentially accelerating the search, into experiments on in vitro evolution. We extend the key instrument of the Holland's theory, the Royal Roads fitness function, to problems of the in vitro evolution (Biological Royal Staircase, BioRS, functions). The specific version of BioRS developed in this publication arises from the realities of experimental evolutionary search for (DNA-) RNA-devices (aptazymes). Our numerical tests showed that for problems with the BioRS functions, simple heuristic algorithms, which turned out to be very effective for preserving BBs in GA, can be very effective in in vitro evolution approaches. We are convinced that such algorithms can be implemented in modern methods of in vitro evolution to achieve significant savings in time and resources and a significant increase in the efficiency of evolutionary search.

Quantification reveals early dynamics in Drosophila maternal gradients.

The Bicoid (Bcd) protein is a primary determinant of early anterior-posterior (AP) axis specification in Drosophila embryogenesis. This morphogen is spatially distributed in an anterior-high gradient, and affects particular AP cell fates in a concentration-dependent manner. The early distribution and dynamics of the bicoid (bcd) mRNA, the source for the Bcd protein gradient, is not well understood, leaving a number of open questions for how Bcd positional information develops and is regulated. Confocal microscope images of whole early embryos, stained for bcd mRNA or the Staufen (Stau) protein involved in its transport, were processed to extract quantitative AP intensity profiles at two depths (apical-under the embryo surface but above the nuclear layer; and basal-below the nuclei). Each profile was quantified by a two- (or three-) exponential equation. The parameters of these equations were used to analyze the early developmental dynamics of bcd. Analysis of 1D profiles was compared with 2D intensity surfaces from the same images. This approach reveals strong early changes in bcd and Stau, which appear to be coordinated. We can unambiguously discriminate three stages in early development using the exponential parameters: pre-blastoderm (1-9 cleavage cycle, cc), syncytial blastoderm (10-13 cc) and cellularization (from 14A cc). Key features which differ in this period are how fast the first exponential (anterior component) of the apical profile drops with distance and whether it is higher or lower than the basal first exponential. We can further discriminate early and late embryos within the pre-blastoderm stage, depending on how quickly the anterior exponential drops. This relates to the posterior-wards spread of bcd in the first hour of development. Both bcd and Stau show several redistributions in the head cytoplasm, quite probably related to nuclear activity: first shifting inwards towards the core plasm, forming either protrusions (early pre-blastoderm) or round aggregations (early nuclear cleavage cycles, cc, 13 and 14), then moving to the embryo surface and spreading posteriorly. These movements are seen both with the 2D surface study and the 1D profile analysis. The continued spreading of bcd can be tracked from the time of nuclear layer formation (later pre-blastoderm) to the later syncytial blastoderm stages by the progressive loss of steepness of the apical anterior exponential (for both bcd and Stau). Finally, at the beginning of cc14 (cellularization stage) we see a distinctive flip from the basal anterior gradient being higher to the apical gradient being higher (for both bcd and Stau). Quantitative analysis reveals substantial (and correlated) bcd and Stau redistributions during early development, supporting that the distribution and dynamics of bcd mRNA are key factors in the formation and maintenance of the Bcd protein morphogenetic gradient. This analysis reveals the complex and dynamic nature of bcd redistribution, particularly in the head cytoplasm. These resemble observations in oogenesis; their role and significance have yet to be clarified. The observed co-localization during redistribution of bcd and Stau may indicate the involvement of active transport.

Quantitative Analysis of the Dynamics of Maternal Gradients in the Early Drosophila Embryo.

Predetermination, formation, and maintenance of the primary morphogenetic gradient (bicoid, bcd) of the early Drosophila embryo involves many interrelated processes. Here we focus on the biological systems analysis of the bcd mRNA redistribution in an early embryo. The results of the quantitative analysis of experimental data, together with the results of their dynamic modeling, substantiate the role of active transport in the redistribution of the bcd mRNA. The role of the nonlinearity of degradation mechanisms in the mRNA pattern robustness is discussed.

Modeling SELEX for regulatory regions using Royal Road and Royal Staircase fitness functions.

The field of evolutionary algorithms (EAs) emerged in the area of computer science due to transfer of ideas from biology and developed independently for several decades, enriched with techniques from probability theory, complexity theory and optimization methods. In this paper, we consider some recent results form the EAs theory transferred back into biology. The well-known biotechnological procedure SELEX (Systematic Evolution of Ligands by EXponential enrichment) is viewed as an experimental implementation of an evolutionary algorithm. Theoretical bounds on EAs runtime are applied to model SELEX search for a regulatory region consisting of promoter and enhancer sequences. A comparison of theoretical bounds to the results of computational simulation indicates some cases where the theoretical bounds give favorable prediction, while simulation requires prohibitive computational resource.

Gene regulatory networks in Drosophila early embryonic development as a model for the study of the temporal identity of neuroblasts.

Genes belonging to the "gap" and "gap-like" family constitute the best-studied gene regulatory networks (GRNs) in Drosophila embryogenesis. Gap genes are a core of two subnetworks controlling embryonic segmentation: (hunchback, hb; Krüppel, Kr; giant, gt; and knirps, kni) and (hb; Kr; pou-domain, pdm; and, probably, castor, cas). Of particular interest is that (hb, Kr, pdm, cas) also specifies the temporal identity of stem cells, neuroblasts, in Drosophila neurogenesis. This GRN controls the sequential differentiation of neuroblasts during the asymmetric cell division. In the last decades, modeling of the patterning of gene ensemble (hb, Kr, gt, kni) in segmentation was in the center of attention. We show that our previously published and extensively studied model at a certain level of external factors is able to reproduce temporal patterns of (hb, Kr, pdm, cas) in neurogenesis with minor evolutionary explicable modifications. This result testifies in favor of a hypothesis that the similarity of two gene ensembles active in segmentation and neurogenesis is a result of co-option of the network architecture in evolution from the common ancestral form. By means of the model dynamical analysis, it is shown that the establishment of the robust patterns in both systems could be explained in terms of the action of attractors in the gap gene dynamical system. We formulate the common principles underlying the robustness of both GRNs in segmentation and neurogenesis due to the similar functional organization of the gene ensembles as having the same evolutionary origin.

Noise model estimation with application to gene expression.

Algorithms for the estimation of noise level and the detection of noise model are proposed. They are applied to gene expression data for Drosophila embryos. The 2D data on gene expression and the extracted 1D profiles are considered. Since the 1D data contain processing errors, an algorithm for separation of these processing errors is constructed to estimate the biological noise level. An approach to discrimination between the additive and multiplicative models is suggested for the 1D and 2D cases. Singular spectrum analysis and its 2D extension are exploited for the pattern extraction. The algorithms are tested on artificial data similar to the real data. Comparison of the results, which are obtained by the 1D and 2D methods, is performed for Krüppel and giant genes.

Two-Exponential Models of Gene Expression Patterns for Noisy Experimental Data.

Spatial pattern formation of the primary anterior-posterior morphogenetic gradient of the transcription factor Bicoid (Bcd) has been studied experimentally and computationally for many years. Bcd specifies positional information for the downstream segmentation genes, affecting the fly body plan. More recently, a number of researchers have focused on the patterning dynamics of the underlying bcd messenger RNA (mRNA) gradient, which is translated into Bcd protein. New, more accurate techniques for visualizing bcd mRNA need to be combined with quantitative signal extraction techniques to reconstruct the bcd mRNA distribution. Here, we present a robust technique for quantifying gradients with a two-exponential model. This approach (1) has natural, biologically relevant parameters and (2) is invariant to linear transformations of the data arising due to variation in experimental conditions (e.g., microscope settings, nonspecific background signal). This allows us to quantify bcd mRNA gradient variability from embryo to embryo (important for studying the robustness of developmental regulatory networks); sort out atypical gradients; and classify embryos to developmental stage by quantitative gradient parameters.

Correction: Cortical movement of Bicoid in early Drosophila embryos is actin- and microtubule-dependent and disagrees with the SDD diffusion model.

[This corrects the article DOI: 10.1371/journal.pone.0185443.].

Relative sensitivity analysis of the predictive properties of sloppy models.

Commonly among the model parameters characterizing complex biological systems are those that do not significantly influence the quality of the fit to experimental data, so-called "sloppy" parameters. The sloppiness can be mathematically expressed through saturating response functions (Hill's, sigmoid) thereby embodying biological mechanisms responsible for the system robustness to external perturbations. However, if a sloppy model is used for the prediction of the system behavior at the altered input (e.g. knock out mutations, natural expression variability), it may demonstrate the poor predictive power due to the ambiguity in the parameter estimates. We introduce a method of the predictive power evaluation under the parameter estimation uncertainty, Relative Sensitivity Analysis. The prediction problem is addressed in the context of gene circuit models describing the dynamics of segmentation gene expression in Drosophila embryo. Gene regulation in these models is introduced by a saturating sigmoid function of the concentrations of the regulatory gene products. We show how our approach can be applied to characterize the essential difference between the sensitivity properties of robust and non-robust solutions and select among the existing solutions those providing the correct system behavior at any reasonable input. In general, the method allows to uncover the sources of incorrect predictions and proposes the way to overcome the estimation uncertainties.

Robustness of expression pattern formation due to dynamic equilibrium in gap gene system of an early Drosophila embryo.

The first manifestation of a segmentation pattern in the early Drosophila development is the formation of expression domains of genes belonging to the gap class. In our previous research the phenomenon of the gap system's robustness, exhibited as the ability to reduce highly variable gene expression in the course of development, was explained as a result of gene cross-regulation. In this paper we formulate the rigorous robustness conditions using the inherent properties of gap gene family. We propose an approach based on the observation that for the formation of a pattern with well-established domain borders it is necessary that there exist a stationary nucleus in which the gene expression level is almost constant in time. The dynamics of gap gene expression is described by a gene circuit model that correctly reproduces the observed principles of the border positioning. We take an advantage of the fact that the system behavior in a stationary nucleus is described by an algebraic equation and hence can be easily handled analytically. This enables us to explicitly characterize the gene cross-regulation properties guaranteeing the system robustness through the spatial behavior of the patterning system along the main embryo axis. In particular, it is proved that if the total regulatory action of all the genes acting in the border formation area changes synchronously with the initial gene gradient, the system will filtrate the initial positioning error and action of highly variable external input factors. We will now show how this approach can be applied as an instrument for the analysis of the robustness mechanism and revealing new biological findings.

Cortical movement of Bicoid in early Drosophila embryos is actin- and microtubule-dependent and disagrees with the SDD diffusion model.

The Bicoid (Bcd) protein gradient in Drosophila serves as a paradigm for gradient formation in textbooks. The SDD model (synthesis, diffusion, degradation) was proposed to explain the formation of the gradient. The SDD model states that the bcd mRNA is located at the anterior pole of the embryo at all times and serves a source for translation of the Bicoid protein, coupled with diffusion and uniform degradation throughout the embryo. Recently, the ARTS model (active RNA transport, synthesis) challenged the SDD model. In this model, the mRNA is transported at the cortex along microtubules to form a mRNA gradient which serves as template for the production of Bcd, hence little Bcd movement is involved. To test the validity of the SDD model, we developed a sensitive assay to monitor the movement of Bcd during early nuclear cycles. We observed that Bcd moved along the cortex and not in a broad front towards the posterior as the SDD model would have predicted. We subjected embryos to hypoxia where the mRNA remained strictly located at the tip at all times, while the protein was allowed to move freely, thus conforming to an ideal experimental setup to test the SDD model. Unexpectedly, Bcd still moved along the cortex. Moreover, cortical Bcd movement was sparse, even under longer hypoxic conditions. Hypoxic embryos treated with drugs compromising microtubule and actin function affected Bcd cortical movement and stability. Vinblastine treatment allowed the simulation of an ideal SDD model whereby the protein moved throughout the embryo in a broad front. In unfertilized embryos, the Bcd protein followed the mRNA which itself was transported into the interior of the embryo utilizing a hitherto undiscovered microtubular network. Our data suggest that the Bcd gradient formation is probably more complex than previously anticipated.

Transcriptional bursting in Drosophila development: Stochastic dynamics of eve stripe 2 expression.

Anterior-posterior (AP) body segmentation of the fruit fly (Drosophila) is first seen in the 7-stripe spatial expression patterns of the pair-rule genes, which regulate downstream genes determining specific segment identities. Regulation of pair-rule expression has been extensively studied for the even-skipped (eve) gene. Recent live imaging, of a reporter for the 2nd eve stripe, has demonstrated the stochastic nature of this process, with 'bursts' in the number of RNA transcripts being made over time. We developed a stochastic model of the spatial and temporal expression of eve stripe 2 (binding by transcriptional activators (Bicoid and Hunchback proteins) and repressors (Giant and Krüppel proteins), transcriptional initiation and termination; with all rate parameters constrained by features of the experimental data) in order to analyze the noisy experimental time series and test hypotheses for how eve transcription is regulated. These include whether eve transcription is simply OFF or ON, with a single ON rate, or whether it proceeds by a more complex mechanism, with multiple ON rates. We find that both mechanisms can produce long (multi-minute) RNA bursts, but that the short-time (minute-to-minute) statistics of the data is indicative of eve being transcribed with at least two distinct ON rates, consistent with data on the joint activation of eve by Bicoid and Hunchback. We also predict distinct statistical signatures for cases in which eve is repressed (e.g. along the edges of the stripe) vs. cases in which activation is reduced (e.g. by mutagenesis of transcription factor binding sites). Fundamental developmental processes such as gene transcription are intrinsically noisy; our approach presents a new way to quantify and analyze time series data during developmental patterning in order to understand regulatory mechanisms and how they propagate noise and impact embryonic robustness.

Sequential construction of a model for modular gene expression control, applied to spatial patterning of the Drosophila gene hunchback.

Gene network simulations are increasingly used to quantify mutual gene regulation in biological tissues. These are generally based on linear interactions between single-entity regulatory and target genes. Biological genes, by contrast, commonly have multiple, partially independent, cis-regulatory modules (CRMs) for regulator binding, and can produce variant transcription and translation products. We present a modeling framework to address some of the gene regulatory dynamics implied by this biological complexity. Spatial patterning of the hunchback (hb) gene in Drosophila development involves control by three CRMs producing two distinct mRNA transcripts. We use this example to develop a differential equations model for transcription which takes into account the cis-regulatory architecture of the gene. Potential regulatory interactions are screened by a genetic algorithms (GAs) approach and compared to biological expression data.

Shaped 3D singular spectrum analysis for quantifying gene expression, with application to the early zebrafish embryo.

Recent progress in microscopy technologies, biological markers, and automated processing methods is making possible the development of gene expression atlases at cellular-level resolution over whole embryos. Raw data on gene expression is usually very noisy. This noise comes from both experimental (technical/methodological) and true biological sources (from stochastic biochemical processes). In addition, the cells or nuclei being imaged are irregularly arranged in 3D space. This makes the processing, extraction, and study of expression signals and intrinsic biological noise a serious challenge for 3D data, requiring new computational approaches. Here, we present a new approach for studying gene expression in nuclei located in a thick layer around a spherical surface. The method includes depth equalization on the sphere, flattening, interpolation to a regular grid, pattern extraction by Shaped 3D singular spectrum analysis (SSA), and interpolation back to original nuclear positions. The approach is demonstrated on several examples of gene expression in the zebrafish egg (a model system in vertebrate development). The method is tested on several different data geometries (e.g., nuclear positions) and different forms of gene expression patterns. Fully 3D datasets for developmental gene expression are becoming increasingly available; we discuss the prospects of applying 3D-SSA to data processing and analysis in this growing field.

Shaped singular spectrum analysis for quantifying gene expression, with application to the early Drosophila embryo.

In recent years, with the development of automated microscopy technologies, the volume and complexity of image data on gene expression have increased tremendously. The only way to analyze quantitatively and comprehensively such biological data is by developing and applying new sophisticated mathematical approaches. Here, we present extensions of 2D singular spectrum analysis (2D-SSA) for application to 2D and 3D datasets of embryo images. These extensions, circular and shaped 2D-SSA, are applied to gene expression in the nuclear layer just under the surface of the Drosophila (fruit fly) embryo. We consider the commonly used cylindrical projection of the ellipsoidal Drosophila embryo. We demonstrate how circular and shaped versions of 2D-SSA help to decompose expression data into identifiable components (such as trend and noise), as well as separating signals from different genes. Detection and improvement of under- and overcorrection in multichannel imaging is addressed, as well as the extraction and analysis of 3D features in 3D gene expression patterns.

Forced evolution in silico by artificial transposons and their genetic operators: The ant navigation problem.

Modern evolutionary computation utilizes heuristic optimizations based upon concepts borrowed from the Darwinian theory of natural selection. Their demonstrated efficacy has reawakened an interest in other aspects of contemporary biology as an inspiration for new algorithms. However, amongst the many excellent candidates for study, contemporary models of biological macroevolution attract special attention. We believe that a vital direction in this field must be algorithms that model the activity of "genomic parasites", such as transposons, in biological evolution. Many evolutionary biologists posit that it is the co-evolution of populations with their genomic parasites that permits the high efficiency of evolutionary searches found in the living world. This publication is our first step in the direction of developing a minimal assortment of algorithms that simulate the role of genomic parasites. Specifically, we started in the domain of genetic algorithms (GA) and selected the Artificial Ant Problem as a test case. This navigation problem is widely known as a classical benchmark test and possesses a large body of literature. We add new objects to the standard toolkit of GA - artificial transposons and a collection of operators that operate on them. We define these artificial transposons as a fragment of an ant's code with properties that cause it to stand apart from the rest. The minimal set of operators for transposons is a transposon mutation operator, and a transposon reproduction operator that causes a transposon to multiply within the population of hosts. An analysis of the population dynamics of transposons within the course of ant evolution showed that transposons are involved in the processes of propagation and selection of blocks of ant navigation programs. During this time, the speed of evolutionary search increases significantly. We concluded that artificial transposons, analogous to real transposons, are truly capable of acting as intelligent mutators that adapt in response to an evolutionary problem in the course of co-evolution with their hosts.

Mid-embryo patterning and precision in Drosophila segmentation: Krüppel dual regulation of hunchback.

In early development, genes are expressed in spatial patterns which later define cellular identities and tissue locations. The mechanisms of such pattern formation have been studied extensively in early Drosophila (fruit fly) embryos. The gap gene hunchback (hb) is one of the earliest genes to be expressed in anterior-posterior (AP) body segmentation. As a transcriptional regulator for a number of downstream genes, the spatial precision of hb expression can have significant effects in the development of the body plan. To investigate the factors contributing to hb precision, we used fine spatial and temporal resolution data to develop a quantitative model for the regulation of hb expression in the mid-embryo. In particular, modelling hb pattern refinement in mid nuclear cleavage cycle 14 (NC14) reveals some of the regulatory contributions of simultaneously-expressed gap genes. Matching the model to recent data from wild-type (WT) embryos and mutants of the gap gene Krüppel (Kr) indicates that a mid-embryo Hb concentration peak important in thoracic development (at parasegment 4, PS4) is regulated in a dual manner by Kr, with low Kr concentration activating hb and high Kr concentration repressing hb. The processes of gene expression (transcription, translation, transport) are intrinsically random. We used stochastic simulations to characterize the noise generated in hb expression. We find that Kr regulation can limit the positional variability of the Hb mid-embryo border. This has been recently corroborated in experimental comparisons of WT and Kr- mutant embryos. Further, Kr regulation can decrease uncertainty in mid-embryo hb expression (i.e. contribute to a smooth Hb boundary) and decrease between-copy transcriptional variability within nuclei. Since many tissue boundaries are first established by interactions between neighbouring gene expression domains, these properties of Hb-Kr dynamics to diminish the effects of intrinsic expression noise may represent a general mechanism contributing to robustness in early development.