In Silico Biology 3, 0010 (2003); ©2003, Bioinformation Systems e.V.  
BGRS 2002

Parametric stability evaluation in computer experiments on the mathematical model of Drosophila control gene subnetwork

Rustem N. Tchuraev and Alexander V. Galimzyanov*




Institute of Biology, Ufa Research Center, Russian Academy of Sciences, Ufa, Russia
Email: galim@anrb.ru
*  corresponding author





Edited by E. Wingender; received September 30, 2002; revised and accepted February 05, 2003; published March 30, 2003


Abstract

Using the method of generalized threshold models, the problem is formulated and solved to evaluate the parametric stability of the model of a gene subnetwork controlling the early ontogenesis of the fruit fly Drosophila melanogaster. Computer experiments have been performed to test the parametric stability of the model. Quantitative evaluations have been obtained for parametric stability of the Drosophila gene subnetwork in nuclei along the embryo's anterior-posterior axis. The results of computer experiments have been compared with the previous research data on "sensitivity" of functioning regimes to random changes of the parameters in the models of prokaryotic and eukaryotic systems, namely the system controlling the -phage development and the subsystem controlling the flower morphogenesis of Arabidopsis thaliana. The obtained results confirm high parametric stability of gene networks that control the development of organisms.

Key words: control gene networks, Drosophila melanogaster, generalised threshold models, computer experiments, parametric stability


Introduction

Nowadays the very rapid accumulation of experimental data reveal the structure of prokaryotic and eukaryotic gene expression control systems as well as the molecular mechanisms of the underlying regulatory protein/gene interactions. Among them gene networks that control ontogenetic processes occupy a particular place. All living objects possess a hierarchically ordered set of control systems, and one control system is judged to be of a higher hierarchical level in relation to the other if the former's functioning results in a change of the other's state [Lyapunov, 1963]. The systems that control ontogenesis are those of a high hierarchical level.

When investigating by mathematical tools the dynamics of gene networks controlling ontogenetic processes it is necessary to evaluate the structural, and particularly parametric, stability of the constructed models. It is evident indeed that every living system possesses the property of stability, in the broad sense, that affords the processes of self-reproduction. Specifically, ontogenesis as a set of intrasystem processes that result in self-reproduction of living systems must be stable both to external influences and internal fluctuations. Ontogenesis stability is achieved to a great extent by the system controlling ontogenesis. At the molecular level the system controlling ontogenesis is realized by an ensemble of interacting intracellular control gene networks. The intracellular gene networks control the activities of its genes and determine properties of the cells and its developmental fate by means of genes' products. Therefore, ontogenesis stability is available particularly through stability of the intracellular control gene network (hereinafter the control gene network, or CGN).

A system is said to be structurally stable if it is stable to fluctuations in parameters and functions [Arnold, 2000]. In the case of control gene networks (CGN) we understand the structural stability as a gene network ability to preserve its inherent normal regimes of functioning (stationary, transitional, periodic) under modifications of the network structure, i. e. by addition/removal of genes or informational links, and also under changes of values for kinetic parameters. The values of such coefficients as unity intensities of gene transcription, mRNA translation, transcript and protein degradation as well as processing and transportation parameters are determined in primary nucleotide sequences of the respective sites. This is true also for the action thresholds of regulatory proteins and their complexes, since they depend mainly on the degree of affinity to the relevant regulatory regions. Thus, fluctuations in the values of kinetic parameters reflect changes in the primary structure of genome DNA molecules. Suppose that in a relatively wide range of values the sensitivity evaluation of functioning regimes in the model to random fluctuations of the parameters is a evaluation of the system's parametric stability. We accent here that the system's parametric stability is an essential component of the structural stability of the system.

The paper presents the results of computer experiments on testing the parametric stability of eukaryotic molecular genetic systems controlling gene expression. The example used is a mathematical model of the gene subnetwork controlling early ontogenesis of the fruit fly Drosophila melanogaster. The model has been constructed on the basis of the method of generalized threshold models [Tchuraev, 1991] in its realization as a computer program [Galimzyanov, 2000].


Methods


Description of the method of generalized threshold models

To solve the problem of predicting the dynamics of a gene network with the known structure we took the method of generalized threshold models (GTM) [Tchuraev, 1991], which allows one to express correctly, in terms of mathematics, some basic regularities in functioning of the system controlling eukaryotic gene expression. Let us next consider the essence of the GTM formalism in accordance with the above-mentioned work.

Taking into account the peculiarity of control processes at the molecular level, this method makes it possible to get both qualitative and quantitative patterns of the dynamics of gene networks. In the molecular system of encoding polymers (DNA, RNA, proteins) and metabolites (m-system) a control subsystem has been singled out that forms controlling variables showing the action of regulatory molecules and a controlled subsystem where, depending on the values of controlling variables controlled variables are formed, i. e. concentrations of DNA, mRNA, proteins and metabolites. The former subsystem is described in terms of discrete mathematics and the latter in terms of the theory of differential equations with particular right parts. Figure 1 presents the microstructure of the genetic block.



Figure 1: Information microstructure of a genetic block. xi (t) is the concentration of the regulatory substance. Di are the discriminators. The discriminator transforms continuous variables xi(t) into discrete binary variables ei (t) in such a way that ei (t) = H [xi (t) - Pij], where H () is the Heaviside function, Pij is the threshold value of concentration, 1, ..., i ..., l ..., n2 are the designations of "microelements" of elementary subautomata of the whole automation Lj, Sj is the designation of a combinative element (combinator of binary signals), TIj(1) is the delay element, j is the control vector.


In the case that the m-system, i. e. the molecular system of encoding polymers and metabolites, contains one gene copy, and the control of the production of the j-th protein is realized at the transcription level, dynamic equations of the following form have been previously obtained:

(1)

where j, j are the columns of n j x 1 dimension; 1j, 1j are the diagonal matrices of nj x nj dimension; T2j is the row of 1 x nj dimension; rj, b2j are the scalars, and the components of the vector-column j denote the current concentration measured by the number of molecules per cell, of the transcripts of the l-th fraction that contains the j-th cistron; j = j(t) is the control vector formed by the logic element and those of delay of the particular block; the diagonal element a(1j)ll of matrix 1j denotes a unity intensity of the transcription beginning from one copy of the l-th promotor; accordingly, component a(2j)l of row T2j is the unity intensity of the translation from the l-th transcript; the diagonal element b(1j)ll of matrix 1j is the unity intensity (degradation factor) of the l-th transcript; b2j is the degradation factor of the j-th polypeptide.

It should be noted that eukaryotic systems controlling gene activity are different from prokaryotic ones:

For a simple case, when regulatory proteins of the i-th form alone affect a regulatory fragment of eukaryotic gene j consisting of kj binding sites, dynamics equations for controlled variables of the eukaryotic genetic block (eukaryotic control gene network) are as follows:

(2)

where mj(0) = mj(0)(t) is the concentration of pre-mRNA molecules containing a j-th cistron; mj = mj(t) is the concentration of "matured" molecules of the j-th RNA fractions related to ribosomes; the parameter is the interval required for processing, transportation and translation of the j-th mRNA fraction; alj is "the specific unity strength" of a single (l-th) regulatory site in the given block; ulj(t) is the control discrete variable as a control vector component; b1j is the degradation coefficient of the j-th transcript; â2j is the translation unity intensity from the j-th transcript; b2j is the degradation coefficient of the j-th polypeptide; rj is the synthesized polypeptide concentration (system output).

The method of generalized threshold models makes it possible to analyze the dynamics of gene networks at both qualitative and quantitative levels: to identify stationary states among any possible regimes of the CGN behaviour, to derive kinetic curves for molecular components (DNA, RNA, proteins) of the control and controlled subsystems. The complexity of mathematical models, as applied to real gene networks, necessitates the development of special accompanying software to study them - GTM formalism is realized in the original computer program package "Analyzer of the Gene Network Dynamics" (AGENDY), its first version being demonstrated at the BGRS'2000 Conference [Galimzyanov, 2000].


Computer program package "Analyzer of the Gene Network Dynamics"

Package modules.   The computational complex AGENDY consists of six principal modules: CEM, MGCS/MGTM, IDV, ODV, PST, and GENEM database. The CEM module (by "create/edit/modify") provides a computer environment for construction, description and modification of CGN models with arbitrary structural (genetic blocks, regulatory interactions) and functional characteristics. The IDV module (by "input data visualization") ensures visual representation of CGN model's structure and parameters and input of the data by a user. The GTM-method for analysing the dynamics of molecular genetic control system (MGCS) is realised in the MGCS/MGTM module. The ODV module (for "output data visualization") helps to organize a combination of forms for representing output data (kinetic curves of mRNA and protein concentrations, diagrams of gene activities) in table and three-dimensional graphic format. The PST block (for "parametric stability testing") realizes schematic computer experiments to test parametric stability of arbitrary fragments of the gene network models under changes of kinetic parameters from an arbitrary subset. GENEM database (by "gene network model") is intended for storing structural and functional characteristics of CGN models, keeping electronic records of computer experiments and provides for the implementation of all the necessary operations in processing the records: their entry, deletion, selection and automatic distribution in the program environment of the element like CGN model.

Data structure.   In the CEM module the gene network model is realized as a bidirectional connected list that allows one to "combine" its structure dynamically out of an arbitrary number of standard-type elements (Genetic Block, Regulatory Element/Operational Site and Regulatory Interaction). In the GENEM module the representation of CGN model is based on interrelated tables in the Paradox format (MGCS, Genes, Sites). Table MGCS gives general information about each model including designation, qualitative description and identification number; the structure of Genes Table agrees with that of the Genetic Block class; fields of the table "Sites" correspond to those of Regulatory Interaction class, and the combination of the records in this table associated with any one genetic block describes the regulatory zone (operational site) of this block. The conception of describing an object like CGN in the MGCS/MGTM module is based on the mathematical model of control molecular genetic systems (MGCS) in GTM formalism. By describing the elements of the model we have developed the program realization abstract classes (MGCS, Genetic Block, Genetic Block Control System, Genetic Block Controlled System, Time Delay Element, etc.) as a combination of fields, methods and properties that reflect the structure and functional coherence of CGN elements.

The GTM method and its program realization AGENDY are an efficient tool in constructing portrayal models of actual molecular systems controlling gene expression in eukaryotic gene networks.


Objects


Qualitative description of Drosophila early ontogenesis and control gene network fragment

A structural plan of Drosophila body and differentiated structures are already determined at the stage of formation of the cellularized blastoderm. At the close of formation of the embryo's cellular structure blastoderm cells contain not the same sets of active genes and turn out to be determined for further differentiation. Selective expression of genes at the early stages of the embryo's development is related to a great extent with transcription regulation of genes of five classes (maternal coordinate, gap, pair rule, segment polarity and homeotic) that form a gene network controlling Drosophila early ontogenesis (Dr-CGN). Dr-CGN gene network in question consists of 22 genes: bicoid (bcd), caudal (cad) and nanos (nos) of maternal coordinate class, giant (gt), hunchback (hb), knirps (kni), Krüppel (Kr) and tailless (tll) of gap class, even-skipped (eve), fushi tarazu (ftz), hairy (h), odd-skipped (odd), odd-paired (opa), paired (prd) and runt (run) of pair-rule class, engrailed (en), wingless (wg) of the segment-polarity class, abdominal-A (abd-A), Antennapedia (Antp), Deformed (Dfd), Sex combs reduced (Scr), Ultrabithorax (Ubx) of the homeotic class. By means of protein products each of these genes is closely involved in the processes of direct and indirect expression regulation of other genes; also it has sites specific for regulatory molecules of protein products of individual genes in the network and plays a key part in the process of formation of cellularized blastoderm and Drosophila body segmentation.

The successive expression of genes in the network begins immediately after the formation of a zygote with maternal mRNA translation of very early genes bcd, nos, cadmat and hbmat localized at both ends of the egg. The newly synthesized proteins are unevenly distributed along the embryo's anterior-posterior axis and specify the spatial patterns of the expression of gap genes. Gap genes are transcribed in the embryo's superimposed compartments and regulate the activity of each other and also of pair-rule genes at the lower level of the cascade. Pair-rule genes are expressed in the embryo's seven alternative bands that can be observed. They encode transcription factors determining the expression of segment polarity genes. Homeotic genes, a group of late genes, specify the type of each segment.


Mathematical model of Drosophila gene subnetwork (Dr-CGN)

Dr-CGN model includes synthetic blocks. G1, G2, ..., G22 are the blocks of protein products synthesis from genes bcd, nos, gt, hb, kni, Kr, tll, eve, ftz, h, odd, opa, prd, run, en(g), wg, abd-A, Antp, Dfd, Scr, Ubx and cad, respectively. Functional bonds among genes by means of their protein products were determined by the known experimental data given in the information resources for Drosophila genome - The Interactive Fly [Brody, 1999], Hox Pro [Spirov et al., 2000]. Figure 2 shows the schematic block of the model.



Figure 2: Scheme for the gene subnetwork controlling Drosophila melanogaster early ontogenesis. Genetic blocks are grouped together in accordance with genes belonging to their respective classes. The negative regulatory bond among genes (transcription repression or translation inhibition) is denoted with a dotted line, and the positive bond (transcription activation with a solid line.


Each synthetic block is described as a genetic block in terms of GTM formalism in accordance with the general scheme [Tchuraev, 1991]. As an example let us take the description of the block for Kni protein synthesis. The block under consideration has the ordinal number 5 in the set of elements of the Dr-CGN model.

Thus, synthetic block G5, i. e. the block of Kni protein synthesis, is controlled by repressors Gt, Hb, Tll [Capovilla et al., 1992; Hulskamp et al., 1990; Pankratz et al., 1989] and activators Bcd, Cad [Rivera-Pomar, 1995; Schulz and Tautz, 1995]. Let us assign ordinal numbers 1 to 5 to regulatory substances Bcd, Gt, Hb, Tll, Cad (from left to right). Then we put three sets of variables: {x15(t), ..., x55(t)} are the concentrations of regulatory substances (input variables), {P15, ..., P55} are the thresholds of protein action upon the regulatory sites of kni gene, {e15(t), ..., e55(t)} are the variables that show either presence or absence of threshold doses of the respective regulatory proteins at time moment t.

Consider the microstructure of block G5 according to the scheme shown in Fig. 1. In this case j = 5, n1 = 5, n2 = 5, l = 1,5. The part of the complex index denoting the genetic block number in the variables and in designations of the elements is omitted for sake of simplicity.

The discriminator Dl transforms input variable xl(t) into binary variable el(t) according to the formula:

(3)

The outputs of discriminators Dl are connected with the inputs of finite determinate automaton 5. automaton 5 consists of five subautomata 1, ..., 5. automaton l describes the interaction of protein l with the regulatory site of the kni gene. Each elementary subautomaton l has a single input channel, along which binary input signal el(t) comes in and a single output channel along which binary output signal l(t) goes out. Thus input (A) and output (B) alphabets in subautomaton l are as follows: A(l) = {0; 1}, B(l) = {0; 1}. A set of states (C) in subautomaton l: C(l) = { q0l, ..., qll}, where l is the refractivity period equalling the average life-time of the complex regulatory protein/DNA site. The functions of transitions and output have the form:

(4)

Functioning of automaton l is described by the following relations:

(5)

where l is the binary variable corresponding to the self-state (activity) of the genetic block by signal el. The self-state of the genetic block is the state with a block when no regulatory substances effect it. The symbol "¯" means logic negation.

automaton 5 is described by a set of discrete finite subautomata l (l = 1,5).

The outputs of automaton 5 are connected with the inputs of combinator S5, a combinational scheme in which the output variables' values at moment t depend only on the input variables at the same moment. In other words a combinator is a trivial automation which depends on the type of regulatory relations in a concrete G-block to produce an integral reaction under the influence of various regulatory molecules exerted on the functioning of the given G-block. Output binary variable 5(t) of combinator S5 is related to input variables with Boolean function F5(t):

5(t) = F5(t) 

A distinct form of the logic function is determined by interaction "logics" among particular regulatory molecules with a definite gene. According to Rivera-Pomar et al., (1995), let us suppose that genetic block Gkni is active only if Bcd and Cad protein concentrations surpass their action thresholds on the block at one and the same time. In this model the logic functions are presented as a perfect disjunctive normal form. For the block under consideration function F5(t) has the form:

(6)

The element of delay TI5 takes into account a delay in synthesis of substance, connected with a limited rate of the RNA polymerase, transcribing a respective gene. Signal 5(t) comes to the input of time delay element TI5 whose output symbol u5(t) is described by the relation of the form:

(7)

where is the self-state (active/non-active) of the block when it is unaffected by regulatory substances. For prokaryotes the genetic block self-state is 1 (active), for eukaryotes it is 0 (non-active). It should be noted that the parameter differs from parameters l in expression (5), since the latter describes the genetic block self-state by a particular signal l.

The value 5 is obtained by the expression:

(8)

where l5 is the length of the gene transcription unit, 5 is the velocity of RNA polymerase movement on genome DNA.

The output variable u5(t) of time delay element TI5 is the control variable for the execution unit and specifies its functioning regime. The execution unit consists of transmitter TR5(1), block PT5 of mRNA processing and transportation and transmitter TR5(2), i. e. the elements that describe the mechanisms of mRNA and protein synthesis and are specified by a set of piecewise linear differential equations (2).

The output of genetic block G5 is variable r5(t) accounting for the concentration of protein Kni (in molecules per cell) at time moment t. The protein concentrations (rj) expressed in molecules per cell can be compared by means of scaling with experimental data (expressed in mol/litre) averaged by the totality of cells. For the Dr-CGN model the concentration measurement in molecules per cell is done under a somewhat simplifying assumption, because at the stage when the expression of gap and pair-rule genes begins cell membranes are just in the process of formation around the nuclei.

Other genetic blocks of the model are described in the same way. No regulation at a post-transcriptional level is taken into account in the present computer version of the model. Since Cad and Hb proteins are translated in the Drosophila egg from both maternal and zygotic mRNA, we formally introduce two additional blocks of Cadmat and Hbmat protein synthesis. Maternal cad and hb-transcripts are included in the model as a complex with regulatory proteins.

The functioning regime of an arbitrary gene network is determined by its structure, a combination of the values of its kinetic parameters and initial data. When modelling Drosophila early ontogenesis we deal with one and the same network, since the embryo's nuclei contain similar genetic material. The values of genetic parameters were selected from the known experimental data [Jackson et al., 2000; Nasiadka and Krause, 1999; Reinitz and Sharp, 1995]. As the initial data, taking the location of nuclei in the embryo into account, we took concentrations of morphogen proteins distributed unevenly along the egg's anterior-posterior axis. Figure 3 shows the concentration gradients of Bcd, Nos, Cadmat and Hbmat proteins.



Figure 3: Distribution of maternal morphogens along the egg's anterior-posterior axis. The gradients of Bcd and Nos protein concentrations have an exponential form along the egg's anterior-posterior axis [Driever and Nüsslein-Volhard, 1988], Cadmat protein concentration has a nearly linear course [Macdonald and Struhl, 1986], Hbmat protein concentration has a rectilinear form in the embryo's first third and an exponential form in the rest of the egg [Tautz, 1988; Štanojevic, 1989].


The model does not take into account the distribution of the products of Dr-CGN genes along the embryo's upper-lower axis, since in the dorso-ventral determination of the embryo's structure a key part is played by dorsal, twist, snail, decapentaplegic and other genes [Sánchez et al., 1997]. Moreover, for the product of the gene bcd, for example, one of the key genes of the network no asymmetrical dorso-ventral distribution was actually found [Driever and Nüsslein-Volhard, 1988]. The nuclei coordinate about the dorsal and ventral cavities is specified with the aid of a control line consisting of two parts and acting as a lateral equator [Reinitz and Sharp, 1995]. The length of the control line (or the egg's length) is taken to equal 100 nuclei.

To investigate the mathematical model Dr-CGN we used the AGENDY software. Calculations of the dynamics in the gene network model were performed for each of a hundred nuclei along the control line. Based on the functioning of the model we obtained sets of kinetic curves describing the dynamics in the model starting with different initial data. Figures 4a-e give examples of such curves. In this case different states of the network can be related to unique expression patterns in the egg's definite regions, i. e. "to project" the kinetics in the model onto the embryo's spatial compartments. For gap genes we have obtained expression model patterns (Fig. 5). Taking into account the assumptions accepted in the model these expression patterns correspond to the real ones.



Figure 4: Kinetic curves for proteins to describe the dynamics in the gene network model in the Drosophila egg's anterior (a), in one of the nuclei of the first parasegment, PS1 (b), in one of the nuclei of the sixth parasegment, PS6 (c), in one of the nuclei of the ninth parasegment, PS9 (d), in the Drosophila egg's posterior part (e).

Figure 5: Expression patterns of gap genes based on functioning of the Dr-CGN model. The obtained patterns of genetic block activities are formed by the 50th minute since the model starts functioning. Dynamics are calculated step-by-step (each step takes a minute). Note that block Gj activity is derived from the value of control variable uj(t) (1 - active, 0 - non-active). Steady-state concentrations of gene network products are achieved by the 120th minute, as in the case of genetic block suppression, for example, synthesized RNA and protein need time to be degraded. Those nuclei are shown in black in which a relevant gene is expressed (or in white in case of no expression).


Computer experiments

For parametric stability we tested a fragment of the Drosophila gene network model consisting of synthetic blocks that correspond to maternal and gap genes (mg-fragment). A set of kinetic parameters with the expression patterns of Dr-CGN genes corresponding to normal early development in Drosophila embryo is taken as a "good" one in view with the assumptions accepted in the construction of the model (expression control patterns) (Fig. 5). Four types of computer experiments were conducted.

In the first experiment the values of kinetic parameters were set at random over wide intervals with the limits taken from physical reasoning: mRNA and protein half-life, 5 to 30 min; RNA-polymerase planting to promoters, 8 to 14 min-1; translation intensity, 0.5 to 2 min-1 (expected mRNA and protein concentrations are 100 to 450 molecules per cell and 1000 to 10000 molecules per cell, respectively).

In the second and third experiments the values of kinetic parameters were selected at random within the limits of l-percent (l = 10%, 20%, 30%, 40%) deviation from the values of parameters in the "good" set when the Dr-CGN model was functioning under the normal regime: Gt protein half-life, 8 min; Hb, 8 min; Kni, 19 min; RNA-polymerase planting to promoters, 10 min-1; translation intensity, 1.5 min-1 (expected mRNA concentrations are 120 molecules per cell, proteins up to 2000 molecules per cell); threshold concentrations of regulatory substances are 300 to 1800 molecules per cell. In the second experiment we varied unit intensities of mRNA and protein transcription and translation as well as coefficients of mRNA and protein degradation. In the third experiment threshold values of regulatory protein concentrations were selected at random within the limits of l-percent deviation from the "good" set. The fourth experiment combined the conditions of the second and third experiments.

When testing the model for parametric stability the dynamics of the gene network in each of the egg's one hundred nuclei (with the account for its location relative to the morphogen concentration gradient) at the selected arbitrary set of parameters was compared with the dynamics of the gene network in the same nucleus with the normal regime on the basis of total genetic blocks in the active state. The stability of the Dr-CGN model to fluctuations in the values of kinetic parameters was evaluated according to the average of normal regimes of gene network functioning in all the egg's nuclei of the total number of arbitrary choices. The parametric stability coefficient CPS was determined by the following expression:

 

where n is the number of nuclei in the control line (n = 100), ci is the percentage of normal regimes of gene network functioning in the i-th nucleus of the total number of arbitrary choices (parametric stability local coefficient).



Results

For each type of the experiments the computer made 500 random choices. The coefficient CPS is 73.2% in the first experiment. Parametric stability coefficients obtained in the second, third and fourth experiments are given below.


Table 1: Stability coefficients of the Dr-CGN model
No. Coefficient CPS (in %) at different l
  l = 10% l = 20% l = 30% l = 40%
2 97.8 93 88.6 83.5
3 76.2 59.8 51 45.3
4 74.2 57 47.5 41.1


Figures 6a-d show the curves for evaluating the parametric stability of the model's mg-fragment relative to the fluctuation in the values of kinetic parameters in the first, second, third and fourth experiments.



Figure 6: The curves for evaluating the model's parametric stability in the experiments: a) E(1, -); b) E(2, 30%); c) E(3, 30%); d) E(4, 10%). Designations: E(k, l) is the type of the experiment, where k is the number of the experiment, l is the deviation percentage; n. n. are the nuclei numbers (embryo's anterior pole is to the left, and posterior pole is to the right); ci is the local index of the model's parametric stability; B1, B2, B3, B4, B5, B6, B7 are the boundaries between model expression patterns of genes gt and tll, tll and gt, gt and hb, hb and Kr, Kr and kni, kni and gt, gt and tll, respectively.


Let us summarize the main experimental results:

  1. During the second, third and fourth experiments, the model's parametric stability decreases with an increase in the deviation percentage of kinetic parameter values from those of the basic set.
  2. Fluctuations in regulatory protein activity thresholds tend to reduce the model's parametric stability much greater as compared with fluctuations in intensities of gap-gene RNA and protein synthesis and degradation.
  3. In all the experiments with Drosophila gene subnetworks, the embryonal nuclei have compartments different from each other by their quantitative (the number of variants with normal functioning regimes) and qualitative (the form of a curve fragment) parametric stability characteristics. In general, the boundaries between compartments correspond to the boundaries between gene expression regions gt and tll, tll and gt, gt and hb, hb and Kr, Kr and kni, kni and gt, gt and tll, respectively. This can be attributed to gene alternative functioning in the subnetwork. Thus, the Drosophila embryo is non-uniform from the standpoint of gene network parametric stability in nuclei along the egg's anterior-posterior axis.
  4. Reasoning from the form of curves for evaluating parametric stability in the third and fourth experiments (both left and right boundaries in one and the same compartment have points with maximum and minimum values of the coefficient ci), it can be supposed that the less the difference between concentration maximum and minimum values of a certain morphogene (gradient local amplitude) at compartment boundaries is, the higher is the model sensitivity to fluctuations in morphogen concentration threshold values.



Discussion

Problems regarding the parametric stability of a particular eukaryotic gene expression control system were also examined in von Dassow et al. (2000). The authors developed and studied a mathematical model for a Drosophila gene network consisting of segment polarity genes. The model was derived from nonlinear ordinary differential equations with account for parameters such as mRNA and protein synthesis intensities and half-life as well as intercellular interactions. Among 240,000 random choices for sets of parameter values they found 1192 sets which allowed the model to reproduce actual gene expression patterns. Such sets were called "solutions". Computer experiments showed that fluctuations in the values of individual parameters (even at 100- to 1000-fold differences) did not influence the process of embryo segmentation along the anterior-posterior axis. Besides the process was shown to be stable to variations of initial data.

In contrast to the scheme reported by G. von Dassow and his co-authors (2000), in our experiments parameter (not random) values for a "good" set and also the limits of intervals for parameter values to be chosen at random were estimated on physical grounds and on the basis of the model's qualitative analysis. Parametric stability was tested with the sets of parameter values lying within the limits of these intervals. Stable occurrence percentage of normal functioning regimes in each successive cycle (with each new random choice) is the criterion for terminating the procedure of parametric stability evaluation, i. e. the number of random choices of kinetic parameter values. The number of random choices (500) in the experiments was chosen because above this number, the deviations of the local coefficients of parametric stability comprise only fractions of a percent, and so there was no need to further increase the number of random choices.

It is necessary here to consider the problem of early ontogenesis modeling with regard to nuclear mitosis divisions and transportation of mRNA and protein molecules among adjacent compartments. In principle, the process of nuclear cleavage can be taken into account by distributing molecules of gene products among daughter nuclei and setting new initial data. However, the law of such a distribution remains unknown. The simplest variant involves halving the number of molecules in nuclei. Since under division the volume of nuclei halves, concentrations of substances (here the number of substance molecules per nucleus volume) vary insignificantly. Thresholds of regulatory protein activity that reflect, in particular, probabilities of binding the transcription factors and DNA regulatory sites do not essentially vary, too. Therefore, it can be thought that a dynamic relationship among current protein concentrations and thresholds of their activity is still preserved under division. That is why nuclear cleavage is not explicitly taken into account in this version of the model. As for the second problem, it should be noted that existing diffusion equations are the most adequate to describe diffusion processes in a homogeneous medium. As far as we know, there is no body of mathematics developed enough to describe the processes of molecular transportation in heterogeneous multi-component molecular systems, specifically in liquid crystalline systems. Since Drosophila egg cytoplasm is a heterogeneous medium, it is difficult to employ the already existing methods. Assuming that molecular distribution of gap-gene products with a particular gradient presents a stabilizing factor, our experimental results can be regarded as an low estimate of parametric stability for a given system.

Next we compare the presented data with the results gained on the models of -phage development control system [Ratner and Tchuraev, 1978] and the subsystem controlling flower morphogenesis of Arabidopsis thaliana [Tchuraev and Galimzyanov, 2001] which were previously constructed on the basis of the GTM method. The computer experimental data of two types for the model of -phage development control system [Ratner and Tchuraev, 1978] count in favour of stability, in the broad sense, characteristic of the system controlling -phage development. It is important that under fluctuations of the parameters over wide limits the system turned out to be capable of functioning in 37 percent of cases (in one of the two regimes - either lysogenic or lytic). This system was found to be not very sensitive to fluctuations in the parameters, i. e. parametrically stable.

The tests for parametric stability of the gene subnetwork controlling Arabidopsis flower morphogenesis showed that there are such sets of kinetic parameters that with one and the same behaviour at the qualitative level the system can show much difference in its quantitative characteristics (2-10-fold distinctions at the level of molecular components concentrations); under identical initial data the system is stable at transitions to stationary states; all four stationary functioning regimes are realized in the system provided that thresholds could be achieved.

The comparison of the experimental results for all three real gene expression control systems under consideration makes it possible to draw the following conclusions for gene networks controlling prokaryotic and eukaryotic ontogenetic processes:

Control gene networks as subnetworks of a single intracellular network contain tens and hundreds of genes, and the whole genomes of eukaryotic organisms consist of some thousands and tens of thousands of genes [Rubin et al., 2000]. Consequently, for special mathematical methods aimed at the investigation of real gene networks one of the efficiency criteria, together with an adequacy for the objects of inquiry, is the possibility to model the dynamics of large systems on their basis. When solving such problems the main limiting factor is the restriction on the amount of required computational resources. Mathematical and program aids used by us - the method of generalized threshold models as an investigative procedure and the accompanying program software "Analyzer of the Gene Network Dynamics" - meet the desired requirements satisfactorily enough.


Acknowledgements

We thank N. A. Kolchanov, V. A. Likhoshvay, I. I. Titov and J. Reinitz for critical discussion of the data given in the paper reviews. We also thank A. P. Maslova for her valuable assistance.


References