Historical Routes for Diversification of Domesticated Chickpea Inferred from Landrace Genomics.

According to archaeological records, chickpea (Cicer arietinum) was first domesticated in the Fertile Crescent about 10,000 years BP. Its subsequent diversification in Middle East, South Asia, Ethiopia, and the Western Mediterranean, however, remains obscure and cannot be resolved using only archeological and historical evidence. Moreover, chickpea has two market types: "desi" and "kabuli," for which the geographic origin is a matter of debate. To decipher chickpea history, we took the genetic data from 421 chickpea landraces unaffected by the green revolution and tested complex historical hypotheses of chickpea migration and admixture on two hierarchical spatial levels: within and between major regions of cultivation. For chickpea migration within regions, we developed popdisp, a Bayesian model of population dispersal from a regional representative center toward the sampling sites that considers geographical proximities between sites. This method confirmed that chickpea spreads within each geographical region along optimal geographical routes rather than by simple diffusion and estimated representative allele frequencies for each region. For chickpea migration between regions, we developed another model, migadmi, that takes allele frequencies of populations and evaluates multiple and nested admixture events. Applying this model to desi populations, we found both Indian and Middle Eastern traces in Ethiopian chickpea, suggesting the presence of a seaway from South Asia to Ethiopia. As for the origin of kabuli chickpeas, we found significant evidence for its origin from Turkey rather than Central Asia.

Modeling the Flowering Activation Motif during Vernalization in Legumes: A Case Study of M. trancatula.

In many plant species, flowering is promoted by the cold treatment or vernalization. The mechanism of vernalization-induced flowering has been extensively studied in Arabidopsis but remains largely unknown in legumes. The orthologs of the FLC gene, a major regulator of vernalization response in Arabidopsis, are absent or non-functional in the vernalization-sensitive legume species. Nevertheless, the legume integrator genes FT and SOC1 are involved in the transition of the vernalization signal to meristem identity genes, including PIM (AP1 ortholog). However, the regulatory contribution of these genes to PIM activation in legumes remains elusive. Here, we presented the theoretical and data-driven analyses of a feed-forward regulatory motif that includes a vernalization-responsive FT gene and several SOC1 genes, which independently activate PIM and thereby mediate floral transition. Our theoretical model showed that the multiple regulatory branches in this regulatory motif facilitated the elimination of no-sense signals and amplified useful signals from the upstream regulator. We further developed and analyzed four data-driven models of PIM activation in Medicago trancatula in vernalized and non-vernalized conditions in wild-type and fta1-1 mutants. The model with FTa1 providing both direct activation and indirect activation via three intermediate activators, SOC1a, SOC1b, and SOC1c, resulted in the most relevant PIM dynamics. In this model, the difference between regulatory inputs of SOC1 genes was nonessential. As a result, in the M. trancatula model, the cumulative action of SOC1a, SOC1b, and SOC1c was favored. Overall, in this study, we first presented the in silico analysis of vernalization-induced flowering in legumes. The considered vernalization network motif can be supplemented with additional regulatory branches as new experimental data become available.

Mechanisms of Vernalization-Induced Flowering in Legumes.

Vernalization is the requirement for exposure to low temperatures to trigger flowering. The best knowledge about the mechanisms of vernalization response has been accumulated for Arabidopsis and cereals. In Arabidopsis thaliana, vernalization involves an epigenetic silencing of the MADS-box gene FLOWERING LOCUS C (FLC), which is a flowering repressor. FLC silencing releases the expression of the main flowering inductor FLOWERING LOCUS T (FT), resulting in a floral transition. Remarkably, no FLC homologues have been identified in the vernalization-responsive legumes, and the mechanisms of cold-mediated transition to flowering in these species remain elusive. Nevertheless, legume FT genes have been shown to retain the function of the main vernalization signal integrators. Unlike Arabidopsis, legumes have three subclades of FT genes, which demonstrate distinct patterns of regulation with respect to environmental cues and tissue specificity. This implies complex mechanisms of vernalization signal propagation in the flowering network, that remain largely elusive. Here, for the first time, we summarize the available information on the genetic basis of cold-induced flowering in legumes with a special focus on the role of FT genes.

Application of Domain- and Genotype-Specific Models to Infer Post-Transcriptional Regulation of Segmentation Gene Expression in Drosophila.

Unlike transcriptional regulation, the post-transcriptional mechanisms underlying zygotic segmentation gene expression in early Drosophila embryo have been insufficiently investigated. Condition-specific post-transcriptional regulation plays an important role in the development of many organisms. Our recent study revealed the domain- and genotype-specific differences between mRNA and the protein expression of Drosophila hb, gt, and eve genes in cleavage cycle 14A. Here, we use this dataset and the dynamic mathematical model to recapitulate protein expression from the corresponding mRNA patterns. The condition-specific nonuniformity in parameter values is further interpreted in terms of possible post-transcriptional modifications. For hb expression in wild-type embryos, our results predict the position-specific differences in protein production. The protein synthesis rate parameter is significantly higher in hb anterior domain compared to the posterior domain. The parameter sets describing Gt protein dynamics in wild-type embryos and Kr mutants are genotype-specific. The spatial discrepancy between gt mRNA and protein posterior expression in Kr mutants is well reproduced by the whole axis model, thus rejecting the involvement of post-transcriptional mechanisms. Our models fail to describe the full dynamics of eve expression, presumably due to its complex shape and the variable time delays between mRNA and protein patterns, which likely require a more complex model. Overall, our modeling approach enables the prediction of regulatory scenarios underlying the condition-specific differences between mRNA and protein expression in early embryo.

Dynamical Modeling of the Core Gene Network Controlling Transition to Flowering in Pisum sativum.

Transition to flowering is an important stage of plant development. Many regulatory modules that control floral transition are conservative across plants. This process is best studied for the model plant Arabidopsis thaliana. The homologues of Arabidopsis genes responsible for the flowering initiation in legumes have been identified, and available data on their expression provide a good basis for gene network modeling. In this study, we developed several dynamical models of a gene network controlling transition to flowering in pea (Pisum sativum) using two different approaches. We used differential equations for modeling a previously proposed gene regulation scheme of floral initiation in pea and tested possible alternative hypothesis about some regulations. As the second approach, we applied neural networks to infer interactions between genes in the network directly from gene expression data. All models were verified on previously published experimental data on the dynamic expression of the main genes in the wild type and in three mutant genotypes. Based on modeling results, we made conclusions about the functionality of the previously proposed interactions in the gene network and about the influence of different growing conditions on the network architecture. It was shown that regulation of the PIM, FTa1, and FTc genes in pea does not correspond to the previously proposed hypotheses. The modeling suggests that short- and long-day growing conditions are characterized by different gene network architectures. Overall, the results obtained can be used to plan new experiments and create more accurate models to study the flowering initiation in pea and, in a broader context, in legumes.

Multi-trait multi-locus SEM model discriminates SNPs of different effects.

BACKGROUND: There is a plethora of methods for genome-wide association studies. However, only a few of them may be classified as multi-trait and multi-locus, i.e. consider the influence of multiple genetic variants to several correlated phenotypes. RESULTS: We propose a multi-trait multi-locus model which employs structural equation modeling (SEM) to describe complex associations between SNPs and traits - multi-trait multi-locus SEM (mtmlSEM). The structure of our model makes it possible to discriminate pleiotropic and single-trait SNPs of direct and indirect effect. We also propose an automatic procedure to construct the model using factor analysis and the maximum likelihood method. For estimating a large number of parameters in the model, we performed Bayesian inference and implemented Gibbs sampling. An important feature of the model is that it correctly copes with non-normally distributed variables, such as some traits and variants. CONCLUSIONS: We applied the model to Vavilov's collection of 404 chickpea (Cicer arietinum L.) accessions with 20-fold cross-validation. We analyzed 16 phenotypic traits which we organized into five groups and found around 230 SNPs associated with traits, 60 of which were of pleiotropic effect. The model demonstrated high accuracy in predicting trait values.

Genomic Analysis of Vavilov's Historic Chickpea Landraces Reveals Footprints of Environmental and Human Selection.

A defining challenge of the 21st century is meeting the nutritional demands of the growing human population, under a scenario of limited land and water resources and under the specter of climate change. The Vavilov seed bank contains numerous landraces collected nearly a hundred years ago, and thus may contain 'genetic gems' with the potential to enhance modern breeding efforts. Here, we analyze 407 landraces, sampled from major historic centers of chickpea cultivation and secondary diversification. Genome-Wide Association Studies (GWAS) conducted on both phenotypic traits and bioclimatic variables at landraces sampling sites as extended phenotypes resulted in 84 GWAS hits associated to various regions. The novel haploblock-based test identified haploblocks enriched for single nucleotide polymorphisms (SNPs) associated with phenotypes and bioclimatic variables. Subsequent bi-clustering of traits sharing enriched haploblocks underscored both non-random distribution of SNPs among several haploblocks and their association with multiple traits. We hypothesize that these clusters of pleiotropic SNPs represent co-adapted genetic complexes to a range of environmental conditions that chickpea experienced during domestication and subsequent geographic radiation. Linking genetic variation to phenotypic data and a wealth of historic information preserved in historic seed banks are the keys for genome-based and environment-informed breeding intensification.

Multi-trait analysis of domestication genes in Cicer arietinum - Cicer reticulatum hybrids with a multidimensional approach: Modeling wide crosses for crop improvement.

Domestication and subsequent breeding have eroded genetic diversity in the modern chickpea crop by ˜100-fold. Corresponding reductions to trait variation create the need, and an opportunity, to identify and harness the genetic capacity of wild species for crop improvement. Here we analyze trait segregation in a series of wild x cultivated hybrid populations to delineate the genetic underpinnings of domestication traits. Two species of wild chickpea, C. reticulatum and C. echinospermum, were crossed with the elite, early flowering C. arietinum cultivar ICCV96029. KASP genotyping of F2 parents with an FT-linked molecular marker enabled selection of 284 F3 families with reduced phenological variation: 255 F3 families of C. arietinum x reticulatum (AR) derived from 17 diverse wild parents and 29 F3 families of C. arietinum x echinospermum (AE) from 3 wild parents. The combined 284 lineages were genotyped using a genotyping-by-sequencing strategy and phenotyped for agronomic traits. 50 QTLs in 11 traits were detected from AR and 35 QTLs in 10 traits from the combined data. Using hierarchical clustering to assign traits to six correlated groups and mixed model based multi-trait mapping, four pleiotropic loci were identified. Bayesian analysis further identified four inter-trait relationships controlling the duration of vegetative growth and seed maturation, for which the underlying pleiotropic genes were mapped. A random forest approach was used to explore the most extreme trait differences between AR and AE progenies, identifying traits most characteristic of wild species origin. Knowledge of the genomic basis of traits that segregate in wild-cultivated hybrid populations will facilitate chickpea improvement by linking genetic and phenotypic variation in a quantitative genetic framework.

Prediction of deleterious mutations in coding regions of mammals with transfer learning.

The genomes of mammals contain thousands of deleterious mutations. It is important to be able to recognize them with high precision. In conservation biology, the small size of fragmented populations results in accumulation of damaging variants. Preserving animals with less damaged genomes could optimize conservation efforts. In breeding of farm animals, trade-offs between farm performance versus general fitness might be better avoided if deleterious mutations are well classified. In humans, the problem of such a precise classification has been successfully solved, in large part due to large databases of disease-causing mutations. However, this kind of information is very limited for other mammals. Here, we propose to better use information available on human mutations to enable classification of damaging mutations in other mammalian species. Specifically, we apply transfer learning-machine learning methods-improving small dataset for solving a focal problem (recognizing damaging mutations in our companion and farm animals) due to the use of much large datasets available for solving a related problem (recognizing damaging mutations in humans). We validate our tools using mouse and dog annotated datasets and obtain significantly better results in companion to the SIFT classifier. Then, we apply them to predict deleterious mutations in cattle genomewide dataset.

Dynamical Modeling of the Core Gene Network Controlling Flowering Suggests Cumulative Activation From the FLOWERING LOCUS T Gene Homologs in Chickpea.

Initiation of flowering moves plants from vegetative to reproductive development. The time when this transition happens (flowering time), an important indicator of productivity, depends on both endogenous and environmental factors. The core genetic regulatory network canalizing the flowering signals to the decision to flower has been studied extensively in the model plant Arabidopsis thaliana and has been shown to preserve its main regulatory blocks in other species. It integrates activation from the FLOWERING LOCUS T (FT) gene or its homologs to the flowering decision expressed as high expression of the meristem identity genes, including AP1. We elaborated a dynamical model of this flowering gene regulatory network and applied it to the previously published expression data from two cultivars of domesticated chickpea (Cicer arietinum), obtained for two photoperiod durations. Due to a large number of free parameters in the model, we used an ensemble approach analyzing the model solutions at many parameter sets that provide equally good fit to data. Testing several alternative hypotheses about regulatory roles of the five FT homologs present in chickpea revealed no preference in segregating individual FT copies as singled-out activators with their own regulatory parameters, thus favoring the hypothesis that the five genes possess similar regulatory properties and provide cumulative activation in the network. The analysis reveals that different levels of activation from AP1 can explain a small difference observed in the expression of the two homologs of the repressor gene TFL1. Finally, the model predicts highly reduced activation between LFY and AP1, thus suggesting that this regulatory block is not conserved in chickpea and needs other mechanisms. Overall, this study provides the first attempt to quantitatively test the flowering time gene network in chickpea based on data-driven modeling.

A Pipeline for Classifying Deleterious Coding Mutations in Agricultural Plants.

The impact of deleterious variation on both plant fitness and crop productivity is not completely understood and is a hot topic of debates. The deleterious mutations in plants have been solely predicted using sequence conservation methods rather than function-based classifiers due to lack of well-annotated mutational datasets in these organisms. Here, we developed a machine learning classifier based on a dataset of deleterious and neutral mutations in Arabidopsis thaliana by extracting 18 informative features that discriminate deleterious mutations from neutral, including 9 novel features not used in previous studies. We examined linear SVM, Gaussian SVM, and Random Forest classifiers, with the latter performing best. Random Forest classifiers exhibited a markedly higher accuracy than the popular PolyPhen-2 tool in the Arabidopsis dataset. Additionally, we tested whether the Random Forest, trained on the Arabidopsis dataset, accurately predicts deleterious mutations in Orýza sativa and Pisum sativum and observed satisfactory levels of performance accuracy (87% and 93%, respectively) higher than obtained by the PolyPhen-2. Application of Transfer learning in classifiers did not improve their performance. To additionally test the performance of the Random Forest classifier across different angiosperm species, we applied it to annotate deleterious mutations in Cicer arietinum and validated them using population frequency data. Overall, we devised a classifier with the potential to improve the annotation of putative functional mutations in QTL and GWAS hit regions, as well as for the evolutionary analysis of proliferation of deleterious mutations during plant domestication; thus optimizing breeding improvement and development of new cultivars.

Analysis of Gene Expression Variance in Schizophrenia Using Structural Equation Modeling.

Schizophrenia (SCZ) is a psychiatric disorder of unknown etiology. There is evidence suggesting that aberrations in neurodevelopment are a significant attribute of schizophrenia pathogenesis and progression. To identify biologically relevant molecular abnormalities affecting neurodevelopment in SCZ we used cultured neural progenitor cells derived from olfactory neuroepithelium (CNON cells). Here, we tested the hypothesis that variance in gene expression differs between individuals from SCZ and control groups. In CNON cells, variance in gene expression was significantly higher in SCZ samples in comparison with control samples. Variance in gene expression was enriched in five molecular pathways: serine biosynthesis, PI3K-Akt, MAPK, neurotrophin and focal adhesion. More than 14% of variance in disease status was explained within the logistic regression model (C-value = 0.70) by predictors accounting for gene expression in 69 genes from these five pathways. Structural equation modeling (SEM) was applied to explore how the structure of these five pathways was altered between SCZ patients and controls. Four out of five pathways showed differences in the estimated relationships among genes: between KRAS and NF1, and KRAS and SOS1 in the MAPK pathway; between PSPH and SHMT2 in serine biosynthesis; between AKT3 and TSC2 in the PI3K-Akt signaling pathway; and between CRK and RAPGEF1 in the focal adhesion pathway. Our analysis provides evidence that variance in gene expression is an important characteristic of SCZ, and SEM is a promising method for uncovering altered relationships between specific genes thus suggesting affected gene regulation associated with the disease. We identified altered gene-gene interactions in pathways enriched for genes with increased variance in expression in SCZ. These pathways and loci were previously implicated in SCZ, providing further support for the hypothesis that gene expression variance plays important role in the etiology of SCZ.

dropEst: pipeline for accurate estimation of molecular counts in droplet-based single-cell RNA-seq experiments.

Recent single-cell RNA-seq protocols based on droplet microfluidics use massively multiplexed barcoding to enable simultaneous measurements of transcriptomes for thousands of individual cells. The increasing complexity of such data creates challenges for subsequent computational processing and troubleshooting of these experiments, with few software options currently available. Here, we describe a flexible pipeline for processing droplet-based transcriptome data that implements barcode corrections, classification of cell quality, and diagnostic information about the droplet libraries. We introduce advanced methods for correcting composition bias and sequencing errors affecting cellular and molecular barcodes to provide more accurate estimates of molecular counts in individual cells.

Translating natural genetic variation to gene expression in a computational model of the Drosophila gap gene regulatory network.

Annotating the genotype-phenotype relationship, and developing a proper quantitative description of the relationship, requires understanding the impact of natural genomic variation on gene expression. We apply a sequence-level model of gap gene expression in the early development of Drosophila to analyze single nucleotide polymorphisms (SNPs) in a panel of natural sequenced D. melanogaster lines. Using a thermodynamic modeling framework, we provide both analytical and computational descriptions of how single-nucleotide variants affect gene expression. The analysis reveals that the sequence variants increase (decrease) gene expression if located within binding sites of repressors (activators). We show that the sign of SNP influence (activation or repression) may change in time and space and elucidate the origin of this change in specific examples. The thermodynamic modeling approach predicts non-local and non-linear effects arising from SNPs, and combinations of SNPs, in individual fly genotypes. Simulation of individual fly genotypes using our model reveals that this non-linearity reduces to almost additive inputs from multiple SNPs. Further, we see signatures of the action of purifying selection in the gap gene regulatory regions. To infer the specific targets of purifying selection, we analyze the patterns of polymorphism in the data at two phenotypic levels: the strengths of binding and expression. We find that combinations of SNPs show evidence of being under selective pressure, while individual SNPs do not. The model predicts that SNPs appear to accumulate in the genotypes of the natural population in a way biased towards small increases in activating action on the expression pattern. Taken together, these results provide a systems-level view of how genetic variation translates to the level of gene regulatory networks via combinatorial SNP effects.

Genomic and phenotypic analysis of Vavilov's historic landraces reveals the impact of environment and genomic islands of agronomic traits.

The Vavilov Institute of Plant Genetic Resources (VIR), in St. Petersburg, Russia, houses a unique genebank, with historical collections of landraces. When they were collected, the geographical distribution and genetic diversity of most crops closely reflected their historical patterns of cultivation established over the preceding millennia. We employed a combination of genomics, computational biology and phenotyping to characterize VIR's 147 chickpea accessions from Turkey and Ethiopia, representing chickpea's center of origin and a major location of secondary diversity. Genotyping by sequencing identified 14,059 segregating polymorphisms and genome-wide association studies revealed 28 GWAS hits in potential candidate genes likely to affect traits of agricultural importance. The proportion of polymorphisms shared among accessions is a strong predictor of phenotypic resemblance, and of environmental similarity between historical sampling sites. We found that 20 out of 28 polymorphisms, associated with multiple traits, including days to maturity, plant phenology, and yield-related traits such as pod number, localized to chromosome 4. We hypothesize that selection and introgression via inadvertent hybridization between more and less advanced morphotypes might have resulted in agricultural improvement genes being aggregated to genomic 'agro islands', and in genotype-to-phenotype relationships resembling widespread pleiotropy.

In silico evolution of the Drosophila gap gene regulatory sequence under elevated mutational pressure.

BACKGROUND: Cis-regulatory sequences are often composed of many low-affinity transcription factor binding sites (TFBSs). Determining the evolutionary and functional importance of regulatory sequence composition is impeded without a detailed knowledge of the genotype-phenotype map. RESULTS: We simulate the evolution of regulatory sequences involved in Drosophila melanogaster embryo segmentation during early development. Natural selection evaluates gene expression dynamics produced by a computational model of the developmental network. We observe a dramatic decrease in the total number of transcription factor binding sites through the course of evolution. Despite a decrease in average sequence binding energies through time, the regulatory sequences tend towards organisations containing increased high affinity transcription factor binding sites. Additionally, the binding energies of separate sequence segments demonstrate ubiquitous mutual correlations through time. Fewer than 10% of initial TFBSs are maintained throughout the entire simulation, deemed 'core' sites. These sites have increased functional importance as assessed under wild-type conditions and their binding energy distributions are highly conserved. Furthermore, TFBSs within close proximity of core sites exhibit increased longevity, reflecting functional regulatory interactions with core sites. CONCLUSION: In response to elevated mutational pressure, evolution tends to sample regulatory sequence organisations with fewer, albeit on average, stronger functional transcription factor binding sites. These organisations are also shaped by the regulatory interactions among core binding sites with sites in their local vicinity.

A new stochastic model for subgenomic hepatitis C virus replication considers drug resistant mutants.

As an RNA virus, hepatitis C virus (HCV) is able to rapidly acquire drug resistance, and for this reason the design of effective anti-HCV drugs is a real challenge. The HCV subgenomic replicon-containing cells are widely used for experimental studies of the HCV genome replication mechanisms, for drug testing in vitro and in studies of HCV drug resistance. The NS3/4A protease is essential for virus replication and, therefore, it is one of the most attractive targets for developing specific antiviral agents against HCV. We have developed a stochastic model of subgenomic HCV replicon replication, in which the emergence and selection of drug resistant mutant viral RNAs in replicon cells is taken into account. Incorporation into the model of key NS3 protease mutations leading to resistance to BILN-2061 (A156T, D168V, R155Q), VX-950 (A156S, A156T, T54A) and SCH 503034 (A156T, A156S, T54A) inhibitors allows us to describe the long term dynamics of the viral RNA suppression for various inhibitor concentrations. We theoretically showed that the observable difference between the viral RNA kinetics for different inhibitor concentrations can be explained by differences in the replication rate and inhibitor sensitivity of the mutant RNAs. The pre-existing mutants of the NS3 protease contribute more significantly to appearance of new resistant mutants during treatment with inhibitors than wild-type replicon. The model can be used to interpret the results of anti-HCV drug testing on replicon systems, as well as to estimate the efficacy of potential drugs and predict optimal schemes of their usage.

Dynamics of miRNA driven feed-forward loop depends upon miRNA action mechanisms.

BACKGROUND: We perform the theoretical analysis of a gene network sub-system, composed of a feed-forward loop, in which the upstream transcription factor regulates the target gene via two parallel pathways: directly, and via interaction with miRNA. RESULTS: As the molecular mechanisms of miRNA action are not clear so far, we elaborate three mathematical models, in which miRNA either represses translation of its target or promotes target mRNA degradation, or is not re-used, but degrades along with target mRNA. We examine the feed-forward loop dynamics quantitatively at the whole time interval of cell cycle. We rigorously proof the uniqueness of solutions to the models and obtain the exact solutions in one of them analytically. CONCLUSIONS: We have shown that different mechanisms of miRNA action lead to a variety of types of dynamical behavior of feed-forward loops. In particular, we found that the ability of feed-forward loop to dampen fluctuations introduced by transcription factor is the model and parameter dependent feature. We also discuss how our results could help a biologist to infer the mechanism of miRNA action.

Mechanisms of developmental robustness.

We present a review of noise buffering mechanisms responsible for developmental robustness. We focus on functions of chaperone Hsp90, miRNA, and cross-regulation of gap genes in Drosophila. The noise buffering mechanisms associated with these functions represent specific examples of the developmental canalization, reducing the phenotypical variability in presence of either genetic or environmental perturbations. We demonstrate that robustness often appears as a function of a network of interacting elements and that the system level approach is needed to understand the mechanisms of noise filtering.

Mechanisms of gap gene expression canalization in the Drosophila blastoderm.

BACKGROUND: Extensive variation in early gap gene expression in the Drosophila blastoderm is reduced over time because of gap gene cross regulation. This phenomenon is a manifestation of canalization, the ability of an organism to produce a consistent phenotype despite variations in genotype or environment. The canalization of gap gene expression can be understood as arising from the actions of attractors in the gap gene dynamical system. RESULTS: In order to better understand the processes of developmental robustness and canalization in the early Drosophila embryo, we investigated the dynamical effects of varying spatial profiles of Bicoid protein concentration on the formation of the expression border of the gap gene hunchback. At several positions on the anterior-posterior axis of the embryo, we analyzed attractors and their basins of attraction in a dynamical model describing expression of four gap genes with the Bicoid concentration profile accounted as a given input in the model equations. This model was tested against a family of Bicoid gradients obtained from individual embryos. These gradients were normalized by two independent methods, which are based on distinct biological hypotheses and provide different magnitudes for Bicoid spatial variability. We showed how the border formation is dictated by the biological initial conditions (the concentration gradient of maternal Hunchback protein) being attracted to specific attracting sets in a local vicinity of the border. Different types of these attracting sets (point attractors or one dimensional attracting manifolds) define several possible mechanisms of border formation. The hunchback border formation is associated with intersection of the spatial gradient of the maternal Hunchback protein and a boundary between the attraction basins of two different point attractors. We demonstrated how the positional variability for hunchback is related to the corresponding variability of the basin boundaries. The observed reduction in variability of the hunchback gene expression can be accounted for by specific geometrical properties of the basin boundaries. CONCLUSION: We clarified the mechanisms of gap gene expression canalization in early Drosophila embryos. These mechanisms were specified in the case of hunchback in well defined terms of the dynamical system theory.