Sequence-based model of gap gene regulatory network
© Kozlov et al.; licensee BioMed Central Ltd. 2014
Published: 19 December 2014
The detailed analysis of transcriptional regulation is crucially important for understanding biological processes. The gap gene network in Drosophila attracts large interest among researches studying mechanisms of transcriptional regulation. It implements the most upstream regulatory layer of the segmentation gene network. The knowledge of molecular mechanisms involved in gap gene regulation is far less complete than that of genetics of the system. Mathematical modeling goes beyond insights gained by genetics and molecular approaches. It allows us to reconstruct wild-type gene expression patterns in silico, infer underlying regulatory mechanism and prove its sufficiency.
We developed a new model that provides a dynamical description of gap gene regulatory systems, using detailed DNA-based information, as well as spatial transcription factor concentration data at varying time points. We showed that this model correctly reproduces gap gene expression patterns in wild type embryos and is able to predict gap expression patterns in Kr mutants and four reporter constructs. We used four-fold cross validation test and fitting to random dataset to validate the model and proof its sufficiency in data description. The identifiability analysis showed that most model parameters are well identifiable. We reconstructed the gap gene network topology and studied the impact of individual transcription factor binding sites on the model output. We measured this impact by calculating the site regulatory weight as a normalized difference between the residual sum of squares error for the set of all annotated sites and for the set with the site of interest excluded.
The reconstructed topology of the gap gene network is in agreement with previous modeling results and data from literature. We showed that 1) the regulatory weights of transcription factor binding sites show very weak correlation with their PWM score; 2) sites with low regulatory weight are important for the model output; 3) functional important sites are not exclusively located in cis-regulatory elements, but are rather dispersed through regulatory region. It is of importance that some of the sites with high functional impact in hb, Kr and kni regulatory regions coincide with strong sites annotated and verified in Dnase I footprint assays.
The detailed analysis of transcriptional regulation is crucial for understanding biological processes, and interest in this problem grows due to new large-scale data acquisition techniques. However despite our expanding knowledge of the biochemistry of gene regulation, we lack a quantitative understanding of this process at a molecular level. We do not understand the mechanism of transcription factor (TF) interactions with adaptor proteins, basal transcriptional machinery and chromatin. We do not know why some cis-regulatory elements (CREs) are modular, while other are scattered over many kilobases of DNA. We cannot effectively predict the aspects of spatiotemporal expression mediated by a particular DNA region and which set of transcription factor binding sites (TFBS) forms a functional CRE.
The segment determination network in Drosophila attracts large interest among researches studying mechanisms of transcriptional regulation. The body of fruit fly consists of repeated morphological units called segments. The borders of segments are demarcated (determined) simultaneously during the blastoderm stage, just before the onset of gastrulation. The segment determination is under control of hierarchical cascade of segmentation genes, most of which are transcriptional regulators. These genes fall into 4 classes. At the bottom of the cascade are the maternal co-ordinate genes bicoid (bcd, one letter code - B) and caudal (cad, one letter code - C). The other groups of genes are gap genes (Kruppel (Kr, one letter code - K), giant (gt, one letter code - G), hunchback (hb, one letter code - H), knirps (kni, one letter code - N), tailless (tll, one letter code - T) and huckebain (hkb, one letter code - J), pair-rule and segment-polarity genes.
There is a large amount of experimental data available about the segment determination system. The gap gene system implements the most upstream regulatory layer of the segmentation gene network. It receives inputs from long-range protein gradients encoded by maternal coordinate genes and establishes discrete territories of gene expression. In this process the gap gene cross-regulation plays important role. The formation of gap gene expression domains is a dynamic process: the domains do not form in one place, but instead in the posterior half of the embryo they shift anteriorly during cleavage cycle 14.
At the molecular level we know the genomic location of many functionally verified CREs, as well as identity and binding affinity of sites for relevant regulating TFs. A wealth of genome scale functional studies provide data on Chip-Seq, RNASeq and DNaseI accessibility measurements. The analysis of these datasets demonstrated that maternal co-ordinate and gap TFs bind to thousands of sites across the Drosophila genome and that the dominant force in cells that modifies the intrinsic DNA specificity of TFs is the inhibition of DNA binding by chromatin . High resolution imaging and image processing techniques provide spatiotemporal data on segmentation gene expression at cellular resolution .
In spite of these efforts we still do not understand the molecular mechanisms involved in gap gene regulation. It is known that the the gap regulatory regions usually contain several CREs driving expression in a precise spatiotemporal pattern and often containing large number of apparent redundant sites for the same TF. Certainly this molecular complexity is important, however the mechanisms underlying it remain elusive.
Mathematical modeling extends the boundaries of genetics and molecular approaches. In general the sufficiency of inferred regulatory mechanism cannot be proven without reconstructing the system ab initio. Currently there is no assay, which accurately reproduces eukaryotic transcription in vitro from well-defined reagents. Mathematical modeling allows us to reconstruct wild type gene expression patterns in silico, to infer underlying regulatory mechanism and prove its sufficiency.
Three major classes of mathematical models have been applied to model regulation in gap gene network: Boolean, differential equation-based and thermodynamic (also termed fractional occupancy) models .
Boolean models represent regulatory relations as logic gates and in the gap gene system they were applied to identify feedback loops which account for topology of gene network at steady-state.
The differential equation based models represent a regulatory network by differential equations, in which a set of molecules such as mRNAs and proteins interact by explicit rules defined in terms of kinetic equations. When applied to the gap gene system these models were able to infer regulatory interactions responsible for formation of the expression domain boundaries, as well as to explain mechanisms for the dynamic anterior shifts of gap domains. It should be noted that the deciphering of the mechanisms of domain motion would be impossible with classic genetic approaches in default of mutants defective for any specific domain shift.
Thermodynamic models rely on simple biophysical descriptions of DNA-protein interactions and statistical physics. They attempt to infer information about gene regulation from the sequences of CREs and the binding affinities of TFs to these elements. This formalism was used to model expression levels in constructs driving reporter gene expression from different gap gene regulatory elements.
It should be noted that all these models have advantages and limitations from the perspective of input data quantity, degree of complexity, and the time interval in which they can model gene expression. Boolean models are suitable to work with large amounts of data produced by genome-wide experiments, but they do not in general consider DNA sequence information. Thermodynamic-based models specifically take into account the features of CREs. However these models provide output for a particular time moment and do not capture the system dynamics. On the contrary differential equation models allow scientists to consider transcriptional regulation over continuous time intervals. The primary limitation of these models is the size of gene network, as the number of parameters rapidly grows with increase of gene number and the problem becomes computationally infeasible. Besides, the differential equation based models usually describe gene interactions in terms of activation/repression and the fine details of transcriptional regulation that thermodynamic-based models offer, are not included.
Evidently, to decipher the molecular mechanisms involved in gap gene regulation we need to understand how genetic information encoded in regulatory elements of these genes is translated into dynamical aspects of gap gene expression. This can be achieved by combining strength of both thermodynamic and differential equation based formalisms. Here we present a new model that provides a dynamical description of gap gene regulatory systems, using detailed DNA-based information, as well as spatial TF concentration data at varying time points. We showed that this model correctly reproduced gap gene expression patterns in wild type embryos and is able to predict gap expression patterns in mutants and reporter constructs.
Results and discussion
Sequence based model of gap gene network
We developed a new model of the gap gene regulatory network which takes as input the affinities of predicted TFBS together with spatial TF concentration data. The output of the model are spatial and temporal patterns of four gap genes hb, Kr, gt, and kni in the form of protein concentration profiles over about one hour of development.
The binding sites for TFs Bcd, Cad, Hb, Gt, Kr, Kni, Tll and Hkb were predicted using position weight matrices (PWMs, see Additional file 1 and Methods). The predicted TFBS affinities were calculated based on the PWM score of the corresponding strongest site as in . The spatial TF concentration data were taken from FlyEx database, which contains data on segmentation gene expression at cycles 13 and 14A of the early embryo development .
Our model consists of two layers. The first layer is a thermodynamic based calculation of the gene activation level. We adopt and modify a method of this calculation presented in . The probability of transcriptional gene activation is assumed to be dependent on the rate of basal transcriptional machinery (BTM) recruitment, which is determined by different probabilities of all possible occupancy states of the regulatory region. Each occupancy state represents a different TF binding configuration on the DNA sequence. As many CREs require mechanisms such as synergy, cooperativity, quenching, and direct repression for proper function [6–10] the model incorporates additional mechanistic features such as short range repression and homotypic cooperativity in transcription factor-DNA binding .
The short range repression, also known as quenching, is a mechanism by which repressors influence activators only if they are bound within a "short range" of the activator binding site [12, 13]. According to this mechanism, a bound repressor cannot interact with the basal complex, but instead leads to a new configuration of the enhancer where its neighborhood in the DNA sequence becomes forbidden to binding by any other TF .
One feature of the model which can be incompatible with the gap gene network is the fact that the type of regulatory action (activation or repression) and its strength for a given TF is the same for all target genes. Previous modeling and experimental results showed that this is not true for gap genes, which may simultaneously exhibit self-activation and repression for other gap genes . Taking this into account, we modified the model to allow different regulatory actions for TFs depending on a target gene, as described in more details in the next section.
Following  we consider that transcriptional output is proportional to the probability of the BTM binding. To model the spatiotemporal dynamics of mRNA and protein synthesis in the early embryo, we write two sets of the reaction-diffusion differential equations [15–17]. We add the delay parameter to account for the average time between events of transcription initiation and corresponding protein synthesis.
We modeled, in one dimension, a region of the blastoderm corresponding to the central midline of the embryo. We consider a time period of cleavage cycles 13 and 14A. Cleavage cycle 14A is about one hour long and is divided into 8 temporal classes (T1-T8) of 6.5 minutes each. The number of nuclei along the A-P axis is doubled when going from c13 to c14. The model was fitted to data on gap protein concentrations from the FlyEx database . Optimization was carried out by differential evolution (DEEP) method [18, 19] to minimize the combined objective function. This function is a sum of the residual sum of squared differences between the model output and data, weighted pattern generating potential and a penalty term, which limits a growth of regulatory weights. The weighted pattern generating potential was proposed in  to account not only for the magnitude of difference between model and data, but also for the direction of change.
In order to further validate that the model is sufficient in data description we constructed a random dataset ("negative control") in which the expression patterns of kni and hb, as well as expression patterns of Kr and gt were shuffled with respect to gene regulatory regions. Consequently, the data the model is fitted to may be considered "nonsense". In this test we hoped that no parameter set could be found making the model output to coincide with "nonsense" data. We noted that a portion of resulting parameter sets has very small affinity constants (K < 10−4) for all TFBS of several TFs, and, hence, these TFs are almost switched off. Evidently, such a situation is not feasible and therefore we removed these parameter sets from further analysis. The mean rms score for the obtained set of parameter vectors was 41.07. The boxplot of the rms scores for biological and negative control data is presented in Figure 2. According to Student's t-test with Welch modification t = 11.26, P − value = 5.101 × 10−15, consequently, the difference in rms mean values is statistically significant.
Gene network topology
In segmentation network a TF can function as both activator and repressor. To account for the possibility of dual regulation we introduced the genetic interconnectivity matrix T ab , which characterizes the action of TF b on gene a. The positive elements of the matrix are statistical weights α A of interaction between bound TF and the BTM, while negative elements correspond to the repressor strength β R . We assume that a bound repressor R acts via the short-range repression mechanism.
Prediction of network topology based on classification of T matrix elements.
The parameter estimates for a representative network.
To understand how reliable our model is we performed the identifiability analysis of the model parameters estimated by fitting to experimental data.
We decide about the sensitivity of the model solution to parameter changes by calculating the confidence intervals for the estimated parameter values (see Methods). This calculation is performed under the assumption that error in data is normally distributed. The error in the gene expression data almost linearly increases with the mean concentration, as it happens for the Poisson distribution. We apply the variance-stabilizing transform to both data and model solution in order to make the error independent of the mean. The parameter estimates found for original objective functional turned out to be also the minimizers for the transformed one.
Estimates and 95% confidence intervals for affinity and cooperativity constants K and ω in a representative network.”
The confidence intervals provide the full information about the parameter estimates only in case of parameter independency, otherwise the intervals are overestimated. Moreover, strong correlation between parameters may lead to their non-identifiability, because a change in one parameter value can be compensated by the appropriate changes of another parameters and, hence, does not significantly influence the solution. It was reported that parameters in the thermodynamic models, for example, affinity constants and cooperativity constants, may be correlated . Because of that we investigated the dependencies between parameters using the collinearity analysis of the sensitivity matrix. This method allows to reveal correlated and hence non-identifiable subsets of parameters.
Two- and three-dimensional subsets of T-matrix elements with collinearity indices higher than 4.
It should be noted that the gene network topology revealed in this work is to a large extent in agreement with experimental evidences , however several disparities exist. In our model Bcd activates Kr in some networks and represses in the others. It was shown that in bcd mutant mothers Kr expression is not reduced but expands anteriorly . This fact leads to proposal that high concentrations of Bcd repress Kr [23, 24], however this effect was later explained by the absence of the anterior gt and hb domains . The activating effect of Bcd on Kr is supported by the fact that Kr expression in reporter constructs is activated by Bcd [26, 27]. The finding that Kr is still expressed in embryos from bcd mutant mothers has been explained by general transcription factor activation  or low levels of Hb [24, 29]. Our analysis does not allow us to make the unambiguous interpretation of the mechanisms of Hkb, Tll and Cad action as these TFs repress and activate target genes in much the same number of networks. It is believed that high concentrations of Cad at the posterior of the embryo activate gap genes. However at about 10 - 15 minutes before gastrulation Cad expression domain refines into a narrow posterior stripe . The posterior hb domain is completely absent in tll mutants [30, 31], that suggests activation of posterior hb by Tll. Some data indicates that Hkb does repress hb, Kr and gt. For example, in hkb mutant embryos the posterior hb domain is unable to retract from the posterior pole . Besides, in embryos mutant for the maternal gene vasa (vas), tll and hkb the Kr domain expands further posterior than in those mutant for vas and tll alone . Finally, in embryos mutant for tll the posterior domain of gt expands less to the posterior pole that in tll;hkb double mutants . An explanation for the model failure to provide unambiguous prediction of the mechanism of Cad, Tll and Hkb action can be found in our analysis of parameter identifiability. This analysis showed that many parameters defining gap gene regulation by Hkb, Cad and Tll are non-identifiable (see Table 3 Table 4 and Figure 3) and therefore we cannot draw any conclusion about these interactions.
Prediction of gap expression in Krmutants and reporter constructs
We use parameters estimated from wild-type expression data to predict in silico gap gene expression patterns in Kr mutants and reporter constructs.
These results convincingly demonstrate that our model is able to correctly predict expression patterns in null mutants and reporter constructs from fits to wild-type data only. This provides an independent proof of model correctness and opens a way for its application for deciphering the mechanisms of transcriptional regulation and gene expression, as will be discussed below.
Contribution of individual TFBS to gap gene expression
Functional genomics studies of animal regulatory regions lead to the hypothesis that transcription factors bind to a majority of genes over a quantitative series of DNA occupancy levels and that the weak regulatory interactions may be of biological significance . Here we use our model to corroborate this idea. Specifically we tried to find the answer on three questions. Are TFBS of small functional impact still important for the model output? Does the correlation between the functional significance of TFBS and its binding affinity exist? Are functional important sites dispersed through regulatory region or predominately located within CREs?
where RSS ref is the residual sum of squares error between the wild type expression data and the model solution for the full set of annotated sites, and RSS mut is the same quantity calculated with the site of interest excluded.
We have calculated the RW for each annotated site and for each gap gene regulatory region. In bioinformatics PWM models are generally used to calculate the BS affinity. We found that the RWs of TFBS show very weak correlation with their PWM score (Spearman's rank correlation coefficient ρ = 0.15, P = 3.5 × 10−6; Pearson correlation coefficient r = 0.17, P = 2.7 × 10 −7). This suggests that the influence of a TFBS on the phenotype is to a great extent explained not by the binding strength per se but by the way the binding sites are involved in the gene regulatory network.
Total number of sites used in the model.
Sites with regulatory the weight wr > 0
The Kr regulatory region contains 10 sites for Gt, Kni, Tll and Cad with RW exceeding the theshold. All these sites fall within different CREs in the region (see Table S2 of Additional File 1). Both Gt and Tll sites within the Kr_730 CRE overlap with annotated DNase I footprint sites (Table 6). The Kr expression in the central domain is produced by Kr_730 or Kr_CD1 elements [26, 27, 45, 46].
We identified 28 sites with RW exceeding the threshold for Bcd, Hb, Kr, Kni and Gt in the gt regulatory region. The (gt_(-1)) CRE drives gt expression in both anterior and posterior domains, while three other CREs reproduce reporter gene expression in the posterior (gt_(-3)) and different anterior domains (gt_(-6), gt_(-10)) [47, 36, 48]. Only 5 of identifed sites are located outside of CREs (see Table S3 of Additional file 1).
The kni regulatory region contains several CREs: kni_(-5) produces anterioventral expresssion, kni_223+64 drives expression in the abdominal region and consists of two discrete sub-elements, kni_(+1) produces expression in both regions. In kni_223+64, the 223-bp sub-element contains one Hb and six Cad TFBS and drives Cad-dependent reporter expression, while the 64-bp sub-element has six binding sites for Bcd and mediates Bcd-dependent expression in the anterior part of the embryo. Interestingly, the anterior expression of the 64-bp element becomes repressed when Hb binds to the 223-bp element . We found 19 sites with RW exceeding threshold for Bcd, Hb, Cad, Kr, Kni and Tll in kni regulatory region (see Table S4 of Additional file 1). Only 6 of these sites are located outside the kni CREs. It is important to note that two sites within kni_(223) sub-element overlap with Cad annotated sites confirmed by DNase I footprint assays (Table 6).
To model the regulatory mechanisms underlying the formation of gap gene expression domains we followed the formalism proposed in  and developed a two-layer model, in which firstly the activation level of each target gene in each embryo nucleus and at each time moment was calculated and at the next step mRNA and protein concentrations for this gene were computed. For calculation of the activation level, we adapted and modified the thermodynamic approach in the form proposed in . We calculate mRNA and protein concentrations by means of differential equations. This innovative approach allowed us to connect the DNA-level information to the system dynamics and thus to overcome a serious limitation of the pure thermodynamic-based models which are static by their nature.
We further modified the method proposed in  by replacing the regulatory parameters α A and β R by the genetic inter-connectivity matrix T ab and introduced the delay parameter τ in our differential equations to account for the average time between events of transcription initiation and corresponding protein synthesis. This makes it possible to translate the elementary regulatory events at the DNA level to the level of gene interactions.
Our modeling approach has clear limitations. The promoter state is calculated by using methods of statistical thermodynamics, while the actual expression products result from this promotor state following the dynamics prescribed by the differential equations. This combination of intrinsically static and dynamical methods in one model is only consistent when there is an evident separation of the timescales of corresponding processes, the equilibration process of TF-DNA binding in the nucleus and the production process of transcribed mRNA and translated protein molecules. Taking into account the complex nature of transcription in eukaryotes, we believe that this assumption is a reasonable approximation for Drosophila genes. One indication for the length of transcription time specific for gap genes is in the fact that the gap gene expression products appear only in late cleavage cycles during the early Drosophila development partially because early cycles are too short for appropriate mRNA maturation . On the other hand, the assumption about equilibrium states of the enhancer binding configurations is also only an approximation. There are clear data showing that such processes as nonspecific binding of TF to DNA and the facilitated diffusion of nonspecifically bound TF to a specific site play their role . Despite the thermodynamic approach proved its efficiency in multiple applications, its proper extension for modeling more dynamic binding configurations seems promising.
The model takes as input the affinities of predicted TFBS together with spatial TF concentration data. The output of the model are spatial and temporal patterns of four gap genes hb, Kr, gt, and kni in the form of protein concentration profiles over about one and a half hour of development.
We used four-fold cross validation test and fitting to random dataset to validate the model and proved its sufficiency in data description. The identifiability analysis showed that most model parameters except of some parameters describing regulation by Tll, Hkb and Cad are well identifiable.
We demonstrated that our model is able to correctly predict expression patterns in Kr null mutants and five different reporter constructs from fits to wild-type data only. This provides an independent proof of model correctness and opens a way for its application for deciphering the mechanisms of transcriptional regulation and gene expression.
We used our model in two ways. Firstly, at the level of gene interactions we reconstructed the gap gene network topology and demonstrated that the basic features of this topology are in agreement with previous modeling results and data from literature .
Secondly, at the DNA level we studied the impact of individual TFBS on the model output. We measured this impact by calculating the site regulatory weight as a normalized difference between the residual sum of squares error for the set of all annotated sites and the set, from which the site of interest was left out. We found that the regulatory weights of TFBS show very weak correlation with their PWM score. This suggests that the influence of a TFBS on the phenotype is to a great extent explained not by the binding strength per se but by the way the binding sites are involved in the gene regulatory network. We also demonstrated that sites with low regulatory weight are important for the model output. This result corroborates the hypothesis about the biological significance of weak regulatory interactions . Our in silico experiments also showed that functional important site are not exclusively located in CREs but are rather dispersed through regulatory region. It is of importance that some of the sites with high functional impact in hb, Kr and kni regulatory regions coincide with strong sites annotated and verified in Dnase I footprint assays.
Transcription factor and gene expression data
We used protein concentrations of transcription factors (referred to as TF in the text) Bcd, Cad, Hb, Gt, Kr, Kni, Tll and Hkb from FlyEx database (http://urchin.spbcas.ru, ) as inputs to the model. This database contains data on segmentation gene expression at the protein level and at discrete time points of cycle 13 and eight time classes (T1-T8) of cycle 14A. To estimate unknown parameters we used the expression patterns of gap genes hb, Kr, kni and gt from the same database. Model predictions were tested using gap gene expression patterns from Kr − embryos obtained from Kr1 loss-of-function allele , as well as reporter constructs driving reporter gene expression from the Kr_CD1, Kr_730, gt_(-3), kni_kd and kni_223+64 CREs (see REDFly database ).
For each of four gap genes hb, Kr, kni and gt we predicted binding sites for Bcd, Cad, Hb, Gt, Kr, Kni, Tll and Hkb in the region spanning 12 Kbp upstream and 6 Kbp downstream of the transcription start site. Transcription factor binding sites were predicted with position weight matrices (PWMs) , which were used to calculate the log-odds score of a site . The PWMs were described in  and can be found at http://www.autosome.ru/iDMMPMM/ (see also the Additional file 1). The PWM thresholds were selected as in .
To model the regulatory mechanisms underlying the formation of gap gene expression domains we adapted the formalism proposed in  and developed a two-layer model, in which firstly the promoter occupancy (activation level) of each target gene in each embryo nucleus and at each time moment was calculated. The second layer of the model is based on differential equations and considers both mRNA and protein synthesis.
where σ is a molecular configuration of the regulatory region for gene a, Q a (σ) is the statistical weight of the interaction between TFs and bound basal transcriptional machinery (BTM), and is the statistical weight of configuration σ for nucleus-time coordinate (i, t), that depends on the concentration of all TFs regulating gene a in nucleus i at time moment t (see  for details).
where v TF is TF concentration, LLR( ) is the log-odd score of a site, computed based on the known PWM of the TF and the background nucleotide distribution, S max is the strongest TFBS and K(S max ) is its binding affinity constant.
where denotes the contribution to statistical weight due to interaction between the TFs bound to sites S i and S j , ω ij is cooperativity constant, d represents the distance between the TF binding sites.
The statistical weight of the interaction between TF and bound BTM Q(S) is the product of the terms corresponding to each bound TF in the configuration. He and coathors  assumed that each TF is either an activator or repressor and proposed and "short-range" mechanism for repression based on the existing experimental work on a few well-characterized or synthetic CREs. A bound activator A interacts with the bound BTM with statistical weight αA >1.
We assume that a bound repressor R does not directly interact with the BTM, but acts via short-range repression mechanism presumably making DNA in its "neighborhood" (defined by a range parameter d R ) inaccessible to binding by any other TF. This configuration with statistical weight scaled by a factor of β R competes with those with the chromatin accessible to activators, thus effectively reducing the occupancy of activators. The parameter β R represents the strength of the repressor and may be interpreted as the equilibrium constant of the reaction that changes the chromatin state from accessible to inaccessible. There is no repression effect when β R is close to 0 while all activator sites are shut down in the neighborhood in the case of β R close to ∞.
In our model we consider that a TF can be activator for one target gene and repressor for another. Therefore we replaced the regulatory parameters α A and β R by the genetic inter-connectivity matrix T ab , which characterizes the action of TF b on gene a. The size of the matrix is N g × (N g + N e ), where N g is the number of gap genes in the model (N g = 4) and N e is the number of external regulatory inputs (bcd, cad, tll and hkb genes, which are not regulated by gap genes, N e = 4). Consequently, β R is translated into negative components of the Tab matrix and positive ones become α A in order to calculate the activation level. This is done for each target gene in each nucleus and for each integration time step.
In the simplest approximation, the target gene expression level is proportional to its activation level . We introduced the delay parameter τ to account for the average time between events of transcription initiation and corresponding protein synthesis, as the model is fitted to gene expression data at the protein level.
where n is the cleavage cycle number, and are maximum synthesis rates,, are the diffusion coefficients, and and are decay rates for protein and mRNA of gene a. The decay rates are related to the mRNA and protein half-lives by . The diffusion term was added to the equation for mRNA to smooth the resulting model output as it was too "spiky" without this term. The parameter is the delay parameter.
We model the dynamics of gap gene expression in cleavage cycles 13 and 14A and in one dimension along the central midline of the embryo. The cycle 14A is divided into eight temporal classes of 6.5 min each. The number of nuclei along the A-P axis is doubled when going from cycle 13 to 14A.
The model was fitted to the protein concentration data for gap genes hb, Kr, gt, and kni from the FlyEx database. Parameter values were optimized by the differential evolution (DEEP) method, described in [18, 19].
The total number of optimized parameters in model (1)-(2) is 68. This includes 32 regulatory parameters T ab , 4 basal machinery constants, 8 binding affinity constants, 8 cooperativity constants, 4 range parameters for short range repression, 4 delay parameters and 8 decay rates. The diffusion constants and synthesis rates were fixed during the optimization.
where g, n and t are gene, nucleus and time point respectively and D is the number of available experimental observations.
were p i and r i are respectively predicted and experimentally observed expression in nucleus i, r max is the maximum level of experimentally observed expression.
This function limits the growth of regulatory parameters, which have very wide ranges.
where the weights were obtained experimentally.
where the Jacobian is the so called sensitivity matrix of size N × m, Fα, m , N−m is an α-quantile of F-distribution with m and N - m degrees of freedom.
The confidence intervals of smaller size correspond to more reliable parameter estimates. In the case when the sign of the parameter estimate provides the most important feature, the estimate is assumed identifiable if the confidence interval is bounded away from zero.
The confidence intervals are overestimated for strongly correlated parameters. Correlation of parameters leads to computational errors since the sensitivity matrix is ill-conditioned.
here λ k is the minimal eigenvalue of the submatrix of the Fisher information matrix. High values of collinearity index mean that the subset of parameters is poorly identifiable because at least two parameters in the subset are interrelated. More details on the methods can be found in .
We are thankful to Prof Sergey Nuzhdin and Prof Alexander M.Samsonov for valuable discussions.The model derivation was supported by the RFBR grant 13-01-00405. Model fitting, validation, adaptation for studies of genetic variability, as well as analysis of model predictions and interpretation of results was done with support of the RSCF grant no.14-14-00302.
Publication of this article has been funded by the BGRS\SB-2014 Organizing Committee and the Russian Science Foundation (RSF).
This article has been published as part of BMC Genomics Volume 15 Supplement 12, 2014: Selected articles from the IX International Conference on the Bioinformatics of Genome Regulation and Structure\Systems Biology (BGRS\SB-2014): Genomics. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/15/S12.
- Kaplan T, Li X-Y, Sabo PJ, Thomas S, Stamatoyannopoulos JA, Biggin MD, Eisen MB: Quantitative models of the mechanisms that control genome-wide patterns of transcription factor binding during early drosophila development. PLoS Genet. 2011, 7 (2): 1001290-10.1371/journal.pgen.1001290.View ArticleGoogle Scholar
- Surkova S, Kosman D, Kozlov K, Manu Myasnikova E, Samsonova A, Spirov A, Vanario-Alonso CE, Samsonova M, Reinitz J: Characterization of the Drosophila segment determination morphome. Developmental Biology. 2008, 313 (2): 844-862. 10.1016/j.ydbio.2007.10.037.PubMedPubMed CentralView ArticleGoogle Scholar
- Ay A, Arnosti DN: Mathematical modeling of gene expression: a guide for the perplexed biologist. Crit Rev Biochem Mol Biol. 2011, 46 (2): 137-151. 10.3109/10409238.2011.556597.PubMedPubMed CentralView ArticleGoogle Scholar
- He X, Samee MAH, Blatti C, Sinha S: Thermodynamics-based models of transcriptional regulation by enhancers: the roles of synergistic activation, cooperative binding and short-range repression. PLoS Comput Biol. 2010, 6 (9): :-View ArticleGoogle Scholar
- Pisarev A, Poustelnikova E, Samsonova M, Reinitz J: FlyEx, the quantitative atlas on segmentation gene expression at cellular resolution. Nucleic Acids Research. 2008, 37: 560-566. [http://nar.oxfordjournals.org/cgi/reprint/gkn717v1.pdf]View ArticleGoogle Scholar
- Simpson-Brose M, Treisman J, Desplan C: Synergy between the Hunchback and Bicoid morphogens is required for anterior patterning in Drosophila. Cell. 1994, 78: 855-865. 10.1016/S0092-8674(94)90622-X.PubMedView ArticleGoogle Scholar
- Arnosti D, Gray S, Barolo S, Zhou J, Levine M: The gap protein Knirps mediates both quenching and direct repression in the Drosophila embryo. The EMBO Journal. 1996, 15: 3659-3666.PubMedPubMed CentralGoogle Scholar
- Gray S, Levine M: Short-range transcriptional repressors mediate both quenching and direct repression within complex loci in Drosophila. Genes and Development. 1996, 10: 700-710. 10.1101/gad.10.6.700.PubMedView ArticleGoogle Scholar
- Hewitt GF, Strunk B, Margulies C, Priputin T, Wang XD, Amey R, Pabst B, Kosman D, Reinitz J, Arnosti DN: Transcriptional repression by the Drosophila Giant protein: Cis element positioning provides an alternative means of interpreting an effector gradient. Development. 1999, 126: 1201-1210.PubMedGoogle Scholar
- Lebrecht D, Foehr M, Smith E, Lopes FJP, Vanario-Alonso CE, Reinitz J, Burz DS, Hanes SD: Bicoid cooperative DNA binding is critical for embryonic patterning in Drosophila. Proceedings of the National Academy of Sciences USA. 2005, 102: 13176-13181. 10.1073/pnas.0506462102.View ArticleGoogle Scholar
- Fakhouri WD, Ay A, Sayal R, Dresch J, Dayringer E, Arnosti DN: Deciphering a transcriptional regulatory code: modeling short-range repression in the Drosophila embryo. Molecular Systems Biology. 2010, 6: 341-PubMedPubMed CentralView ArticleGoogle Scholar
- Kulkarni MM, Arnosti DN: cis-Regulatory logic of short-range transcriptional repression in Drosophila melanogaster. Molecular and Cellular Biology. 2005, 25: 3411-3420. 10.1128/MCB.25.9.3411-3420.2005.PubMedPubMed CentralView ArticleGoogle Scholar
- Arnosti DN, Kulkarni MM: Transcriptional enhancers: Intelligent enhanceosomes or flexible billboards?. Journal of Cellular Biochemistry. 2005, 94 (5): 890-898. 10.1002/jcb.20352.PubMedView ArticleGoogle Scholar
- Jaeger J: The gap gene network. Cellular and Molecular Life Sciences. 2011, 68: 243-274. 10.1007/s00018-010-0536-y.PubMedPubMed CentralView ArticleGoogle Scholar
- Reinitz J, Sharp DH: Mechanism of eve stripe formation. Mechanisms of Development. 1995, 49: 133-158. 10.1016/0925-4773(94)00310-J.PubMedView ArticleGoogle Scholar
- Jaeger J, Surkova S, Blagov M, Janssens H, Kosman D, Kozlov KN, Manu Myasnikova E, Vanario-Alonso CE, Samsonova M, Sharp DH, Reinitz J: Dynamic control of positional information in the early Drosophila embryo. Nature. 2004, 430: 368-371. 10.1038/nature02678.PubMedView ArticleGoogle Scholar
- Kozlov K, Surkova S, Myasnikova E, Reinitz J, Samsonova M: Modeling of gap gene expression in Drosophila Kruppel mutants. PLoS Comput Biol. 2012, 8 (8): 1002635-10.1371/journal.pcbi.1002635.View ArticleGoogle Scholar
- Kozlov K, Samsonov A: Deep - differential evolution entirely parallel method for gene regulatory networks. Journal of Supercomputing. 2011, 57: 172-178. 10.1007/s11227-010-0390-6.PubMedPubMed CentralView ArticleGoogle Scholar
- Kozlov K, Ivanisenko N, Ivanisenko V, Kolchanov N, Samsonova M, Samsonov AM: Enhanced Differential Evolution Entirely Parallel Method for Biomedical Applications. LNCS 7979 V. Malyshkin (Ed.): PaCT 2013. 2013, 409-416.Google Scholar
- Samee MAH, Sinha S: Evaluating thermodynamic models of enhancer activity on cellular resolution gene expression data. Methods. 2013, 62: 79-90. 10.1016/j.ymeth.2013.03.005.PubMedView ArticleGoogle Scholar
- Sachs L: Applied Statistics. 1982, Springer, USView ArticleGoogle Scholar
- Dresch JM, Richards M, Ay A: Thermodynamic modeling of transcription: sensitivity analysis differentiates biological mechanism from mathematical model-induced effects. BMC Systems Biology. 2010, 4 (142): 2-11.Google Scholar
- Gaul U, Jackle H: Pole region-dependent repression of the Drosophila gap gene Krüppel by maternal gene products. Cell. 1987, 51: 549-555. 10.1016/0092-8674(87)90124-3.PubMedView ArticleGoogle Scholar
- Hülskamp M, Pfeifle C, Tautz D: A morphogenetic gradient of Hunchback protein organizes the expression of the gap genes Krüppel and knirps in the early Drosophila embryo. Nature. 1990, 346: 577-580. 10.1038/346577a0.PubMedView ArticleGoogle Scholar
- Kraut R, Levine M: Mutually repressive interactions between the gap genes giant and Krüppel define middle body regions of the Drosophila embryo. Development. 1991, 111: 611-621.PubMedGoogle Scholar
- Hoch M, Schröder C, Seifert E, Jäckle H: Cis-acting control elements for Krüppel expression in the Drosophila embryo. The EMBO Journal. 1990, 9: 2587-2595.PubMedPubMed CentralGoogle Scholar
- Hoch M, Seifert E, Jäckle H: Gene expression mediated by cis-acting sequences of the Krüppel gene in response to the Drosophila morphogens Bicoid and Hunchback. The EMBO Journal. 1991, 10: 2267-2278.PubMedPubMed CentralGoogle Scholar
- Kerrigan LA, Croston GE, Lira LM, Kadonaga JT: Sequence-specific transcriptional antirepression of the Drosophila Krüppel gene by the GAGA factor. The Journal of Biological Chemistry. 1991, 266: 574-582.PubMedGoogle Scholar
- Schulz C, Tautz D: Autonomous concentration-dependent activation and repression of Krüppel by hunchback in the Drosophila embryo. Development. 1994, 120: 3043-3049.PubMedGoogle Scholar
- Casanova J: Pattern formation under the control of the terminal system in the Drosophila embryo. Development. 1990, 110: 621-628.PubMedGoogle Scholar
- Brönner G, Jäckle H: Control and function of terminal gap gene activity in the posterior pole region of the Drosophila embryo. Mechanisms of Development. 1991, 35: 205-211. 10.1016/0925-4773(91)90019-3.PubMedView ArticleGoogle Scholar
- Rivera-Pomar R, Niessing D, Schmidt-Ott U, Gehring WJ, Jackle H: RNA binding and translational suppression by Bicoid. Nature. 1996, 379: 746-749. 10.1038/379746a0.PubMedView ArticleGoogle Scholar
- Furriols M, Casanova J: In and out of torso RTK signaling. The EMBO Journal. 2003, 22: 1947-1952. 10.1093/emboj/cdg224.PubMedPubMed CentralView ArticleGoogle Scholar
- Dequeant M-L, Pourquie O: Segmental patterning of the vertebrate embryonic axis. Nat Rev Genet. 2008, 9: 370-382. 10.1038/nrg2320.PubMedView ArticleGoogle Scholar
- Gallo SM, Gerrard DT, Miner D, Simich M, Des Soye B, Bergman CM, Halfon MS: Redfly v3.0: Toward a comprehensive database of transcriptional regulatory elements in Drosophila. Nucleic Acids Research. 2010, 10: 1-6.Google Scholar
- Schroeder MD, Pearce M, Fak J, Fan H-Q, Unnerstall U, Emberly E, Rajewsky N, Siggia ED, Gaul U: Transcriptional control in the segmentation gene network of Drosophila. PLoS Biology. 2004, 2: 271-10.1371/journal.pbio.0020271.View ArticleGoogle Scholar
- Rivera-Pomar R, Lu X, Perrimon N, Taubert H, Jackle H: Activation of posterior gap gene expression in the Drosophila blastoderm. Nature. 1995, 376: 253-256. 10.1038/376253a0.PubMedView ArticleGoogle Scholar
- Biggin MD: Animal transcription networks as highly connected, quantitative continua. Developmental cell. 2011, 21 (4): 611-626. 10.1016/j.devcel.2011.09.008.PubMedView ArticleGoogle Scholar
- Driever W, Nusslein-Volhard C: The Bicoid protein is a positive regulator of hunchback transcription in the early Drosophila embryo. Nature. 1989, 337: 138-143. 10.1038/337138a0.PubMedView ArticleGoogle Scholar
- Struhl G, Struhl K, Macdonald PM: The gradient morphogen Bicoid is a concentration-dependent transcriptional activator. Cell. 1989, 57: 1259-1273. 10.1016/0092-8674(89)90062-7.PubMedView ArticleGoogle Scholar
- Driever W, Thoma G, Nusslein-Volhard C: Determination of spatial domains of zygotic gene expression in the Drosophila embryo by the affinity of binding sites for the Bicoid morphogen. Nature. 1989, 340: 363-367. 10.1038/340363a0.PubMedView ArticleGoogle Scholar
- Treisman J, Desplan C: The products of the Drosophila gap genes hunchback and Krüppel bind to the hunchback promoters. Nature. 1989, 341: 335-337. 10.1038/341335a0.PubMedView ArticleGoogle Scholar
- Margolis JS, Borowsky ML, Steingrimsson E, Shim CW, Lengyel JA, Posakony JW: Posterior stripe expression of hunchback is driven from two promoters by a common enhancer element. Development. 1995, 121: 3067-3077.PubMedGoogle Scholar
- Lukowitz W, Schröder C, Glaser G, Hulskamp M, Tautz D: Regulatory and coding regions of the segmentation gene hunchback are functionally conserved between Drosophila virilis and Drosophila melanogaster. Mechanisms of Development. 1994, 45: 105-115. 10.1016/0925-4773(94)90024-8.PubMedView ArticleGoogle Scholar
- Hoch M, Gerwin N, Taubert H, Jackie H: Competition for overlapping sites in the regulatory region of the Drosophila gene Krüppel. Science. 1992, 256: 94-97. 10.1126/science.1348871.PubMedView ArticleGoogle Scholar
- Capovilla M, Eldon ED, Pirrotta V: The giant gene of Drosophila encodes a b-ZIP DNA-binding protein that regulates the expression of other segmentation gap genes. Development. 1992, 114: 99-112.PubMedGoogle Scholar
- Berman BP, Nibu Y, Pfeiffer BD, Tomancak P, Celniker SE, Levine M, Rubin GM, Eisen MB: Exploiting transcription factor binding site clustering to identify cis-regulatory modules involved in pattern formation in the Drosophila genome. Proceedings of the National Academy of Sciences USA. 2002, 99: 757-762. 10.1073/pnas.231608898.View ArticleGoogle Scholar
- Segal E, Raveh-Sadka T, Schroeder M, Unnerstall U, Gaul U: Predicting expression patterns from regulatory sequence in Drosophila segmentation. Nature. 2008, 451: 535-540. 10.1038/nature06496.PubMedView ArticleGoogle Scholar
- Dresch JM, Thompson MA, Arnosti DN, Chiu C: Two-Layer Mathematical Modeling of Gene Expression: Incorporating DNA-Level Information and System Dynamics. SIAM J APPL M ATH. 2013, 73 (2): 804-826. 10.1137/120887588.View ArticleGoogle Scholar
- Sharon E, van Dijk D, Kalma Y, Keren L, Manor O, Yakhini Z, Segal E: Probing the effect of promoters on noise in gene expression using thousands of designed sequences. Genome research. 2014, 168773-113.Google Scholar
- Kulakovskiy IV, Makeev VJ: Discovery of dna motifs recognized by transcription factors through integration of different experimental sources. Biophysics. 2009, 54 (6): 667-674. 10.1134/S0006350909060013.View ArticleGoogle Scholar
- Lifanov AP, Makeev VJ, Nazina AG, Papatsenko DA: Homotypic regulatory clusters in drosophila. Genome Research. 2003, 13 (4): 579-588. 10.1101/gr.668403. [http://genome.cshlp.org/content/13/4/579.full.pdf+html]PubMedPubMed CentralView ArticleGoogle Scholar
- Kulakovskiy IV, Boeva VA, Favorov AV, Makeev VJ: Deep and wide digging for binding motifs in ChIP-Seq data. Bioinformatics. 2010, 26 (20): 2622-3. 10.1093/bioinformatics/btq488.PubMedView ArticleGoogle Scholar
- Kulakovskiy IV, Favorov AV, Makeev VJ: Motif discovery and motif finding from genome-mapped DNase footprint data. Bioinformatics. 2009, 25 (18): 2318-25. 10.1093/bioinformatics/btp434.PubMedView ArticleGoogle Scholar
- Struffi P, Corado M, Kaplan L, Yu D, Rushlow C, Small S: Combinatorial activation and concentration-dependent repression of the Drosophila even skipped stripe 3+7 enhancer. Development. 2011, 138: , 4291-4299. 10.1242/dev.065987.PubMedPubMed CentralView ArticleGoogle Scholar
- Ashyraliyev M, Jaeger J, Blom JG: On parameter estimation and determinability for drosophila gap gene circuits. BMC Systems Biology. 2008, 2 (83):Google Scholar
- Bates DM, Watts DG: Nonlinear Regression Analysis and Its Applications. 1988, JOHN WILEY & SONS, INC., New York, NYView ArticleGoogle Scholar
- Myasnikova E, Kozlov KN: Statistical method for estimation of the predictive power of a gene circuit model. J Bioinform Comput Biol. 2014, 12 (2): 1441002-10.1142/S0219720014410029.PubMedView ArticleGoogle Scholar
- Brun R, Reichert P, Kunsch HR: Practical identifiability analysis of large environmental simulation models. Water Resour Res. 2001, 37 (4): 1015-1030. 10.1029/2000WR900350.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.