1 Introduction
During the recent years, the rapid development of the bio-arrays techniques [1] based on isotopic or fluorescent activity of hybridised DNA chips allowed the biologist to give to a grey level peak the signification of an expression rate for the genes studied in the bio-array. If we repeat this acquisition at different times of the cell cycle for different cells of a same tissue, we can calculate correlations between the genes expression rate and hence we are able to make explicit a matrix W called the intergenic interaction matrix, representing the repression and induction influences a gene can exert on other genes.
1.1 The raw data from the bio-array imaging
The first encountered problem with the bio-array image is the noise and we have to low-pass it in order to suppress the high-frequency noise (see Fig. 1). The result of this pre-treatment is a better separation of the isotopic activity peaks, allowing a watershed separation and contouring [2,3]. We can also apply a more accurate segmentation and contouring method called potential-Hamiltonian or ‘Gaussian-stamping’ method. Let us remark that the peaks are about Gaussian, with a weak kurtosis and skewness, allowing in particular the respect of the conservation ‘law’: 2/3 of the activity peak are concentrated into the set of points (x,y) where the Gaussian curvature C(x,y) vanishes, i.e. inside the maximum gradient line, called the characteristic line of the peak. Then it is possible to neglect the part of the peak outside the projection of this line, whose equation is C(x,y)=∂2g/∂x2∂2g/∂y2−(∂2g/∂x∂y)2=0, g(x,y) denoting the grey level function at the pixel of coordinates (x,y).

(a) Raw data; (b) low-pass filtering; (c) watershed segmenting; (d) contouring.
We are thus led to consider the new height function C(x,y) instead of the function g(x,y) and its level line C(x,y)=0. A new algorithm has been proposed in [1] to obtain automatically the characteristic line and, after, by integrating inside the projection of this line on the (x,y) plane and multiplying by 3/2 the obtained result, to standardize the estimated activity in terms of a bio-image with small squares symbolizing in grey levels the degree of hybridisation of the cDNA's expressing the regulation of a glioma tissue (Fig. 2). From such bio-array images acquired in cells of the same tissue at different times of the cell cycle, we can study the interactions between genes by estimating an interaction matrix.

‘Gaussian-stamping’ contours (left) and standardized bio-array image (black rectangle right corresponds to the left image).
2 Some rigorous results about the network attractors
2.1 The interaction matrix
The interaction matrix W is similar to the synaptic weight matrix, which rules the relationships between neurons in a neural network. The general coefficient wik of such an interaction matrix W is equal to +1 (resp −1, 0) if the gene Gk activates (resp inhibits, does not influence) the gene Gi, the state xi of the gene Gi being equal to +1 (resp −1), if it is (resp is not) expressed. In the case of small regulatory genetic systems (called operons), the knowledge of such a matrix W permits to make explicit all possible stationary behaviours of the organisms having the corresponding genome. The change of state of gene Gi between t and t+1 obeys a threshold rule: xi(t+1)=H(∑k=1,nwikxk(t)−bi) or

The lactose operon exhibited by Jacques Monod.
If G is the interaction graph associated to W, then we call connected component C of G each set of genes C such that there is a path between every pair of genes of C along a sequence of arcs of G. A Garden of Eden is a gene receiving no arc, but influencing at least one other gene. A regulon is a connected component of G having exactly one positive (auto-catalysis) and one negative loop, these loops sharing the auto-catalysed node. In the lactose operon (see its G in Fig. 3),
2.2 Some definitions and notations
In the following, we give some definitions about the rigorous mathematical description of the discrete Boolean networks used to describe the genes interaction dynamics, and their associated graph G and matrix W. Then we will present some theoretical results with rigorous proofs only for the first and last results in order to show what kind of reasoning we have to perform and we will refer to [13–16] for complete demonstrations. Let us consider a graph G=(V, E), where V={1,…,n} is the set of nodes and E is the set of arcs. Let
2.3 Relations between positive and negative cycles, and fixed points
In the sequel, we will assume that the graph G is connected, since otherwise one can apply the results to each of connected components of G. In addition, we will suppose with no loss of generality that |Γ−(i)|>0, for all i∈V. Since otherwise, if there exists a node i∈V such that Γ−(i) is empty, then we can assume that the arc (i,i) exists in E; in this way the dynamics of both networks are the same. It evidently follows from this property that there exists at least one circuit C in G (possibly a circuit of the form (i,i)). Finally, we suppose that the graph G and the matrix W have a quasi-minimal structure, that is to say, all (j,i), such as i≠j, belong to E (or equivalently wij≠0, if i≠j), if there exists
The following property will be very useful in the following for characterizing a cycle.
Proposition 1
A cycle C is positive if and only if there exists a vector
Proof
Let C be a positive cycle and i(0) a fixed node belonging to C. Let us enumerate the nodes belonging to C by i(0),i(1),…,i(j), such that
Let C be now a negative cycle, and let us suppose that equation (1) is true, then ∏(j,i)∈Csign(wij)=∏jxj∏ixi=(∏jxj)2, but sign(C)=∏(j,i)∈Csign(wij)<0, which is contradictory.□
Theorem 1
Given N, if all cycles of the incidence graph G are positive, then there exists a vector
Remark
There are two remarkable fixed points having by construction a non-frustration property, that is on each cycle the sign changes of xis are identical to the sign changes of the arcs. For other possible fixed points, there is at least one cycle for which sign changes are frustrated.
Theorem 2
If all circuits of the incidence graph G are negative, then N has no fixed points.
2.4 Minimal regulatory networks
The previous results allow us to characterize some minimal regulatory networks. The following propositions constitute examples of minimal regulatory networks. They solve in part the inverse problem consisting in the description of W only from the knowledge of a phenotypic x observed from bio-array images.
Proposition 2
Let N having n nodes and n connections, a necessary and sufficient condition of existence of a fixed point
Proposition 3
Given a state vector
- (1) wik=αikxixk, where αik⩾0 and, for all i, there exists k(i) such that αik(i)≠0 and
- (2) −|αik(i)|<bi⩽|αik(i)|.
Proposition 4
Let N with n nodes and n+1 connections, a necessary and sufficient condition for existence of an attractor of all points parallel iterated, is a negative circuit and a positive circuit intersecting.
2.5 Fixed points bounds in regulatory networks
Given N, let C be a positive circuit of N, then by Proposition 1, there exists
Lemma 1
Given N and
Theorem 3
If m is the total number of positive circuits of N, then the number of fixed points of N is⩽2m, and this upper bound is reached if and only if for all circuits C of N there does not exist an arc (k,i) in CC ending in C (there is no garden of Eden k pending to C).
Remark
We have to notice that the condition concerns the number of circuits and not of cycles, these last being in general very more numerous.
2.6 Asymptotic mean value for the number of fixed configurations (fixed vectors or limit cycles of vectors) in the case K
Let us consider now a network N having n nodes and 2n connections such as
Lemma 2
For any graph G having m non-oriented edges, the mean number of oriented edges we can define on G from the non-oriented configuration is equal to 4m/3.
Proof
Let us note 〈o〉 the mean number of oriented edges we can construct from a configuration of m non-oriented edges; then, if exactly k from the m non-oriented edges are decomposed into two oriented opposite connections, we have Cmm−k2m−k different ways to dispatch the not double connections into the (m−k) other non-oriented edges; hence we can write:
Theorem 4
If the network N has n nodes and
Proof
Following [19], if the connections of N are random, and if the mean number c of non-oriented edges per node is equal to 3/2, then the random variables Xi equal to the number of disjoint cycles of length i of N are independent and Poissonnian with parameter λ(i)=2i−1/i, if n is sufficiently large. From Lemma 2, we are just in this case, because we have 2n=4m/3 connections, hence m=3n/2 and c=m/n=3/2. Then we have, for the mean number 〈f〉 of fixed configurations of N:
We will now evaluate the expectation A(σ). Each disjoint positive circuit bringing two fixed points (Theorem 3 above), an isolated positive non-circuit cycle bringing also two fixed points and an isolated negative circuit bringing one limit cycle, we can first calculate A(0,σ), the expected number of fixed configurations in the case where we have only disjoint positive circuits, the rest of the nodes depending on these circuits (and hence their states being fixed by the states of the circuit): A(0,σ)=B(0,σ)/D(σ), where B(0,σ)=2s (number of fixed points of s=∑i=1ns(i) disjoint positive circuits, from Theorem 3 above)×2k (number of different signs for each of the k=∑i=1nis(i)⩽n nodes involved in the s circuits)×[2s (number of different directions – left or right – for each of the s circuits)/2s (reduction factor for having only positive circuits)]×N(σ),D(σ)=2k (number of different directions for each of the k connections)×2k (number of different signs for each of the k connections)×N(σ), where N(σ) is the number of choices for the s disjoint cycles: N(σ)=Cks(1)1,…,s(n)n∏i=1n(i−1)!s(i)/2s. N(σ) just equals the number of choices of k nodes in s(1) subsets of size 1,…,s(n) of size n times the number of choices of different loops (without multiple points) connecting the vertices inside each of these subsets.
In the same way, we can calculate A(1,σ) (resp A(j,σ)) the expected number of attractors of N in the case where we have among the s disjoint cycles 1 (resp j) isolated positive non-circuit cycles (bringing two fixed vectors) or isolated negative circuits (bringing 1 fixed cycle).
We have also:
Finally, more generally, we have A(j,σ)=B(j,σ)/D(σ), where
By summing the A(j,σ)s and after the A(σ)Πσs, 〈f〉 is clearly of the order of n1/2:
Then we have:
Remark
3 Examples of genetic regulation networks
3.1 The flowering regulatory network of Arabidopsis thaliana
If we consider the interaction graph of the flowering regulatory network of Arabidopsis thaliana (Fig. 4, right) [10], then we can easily define from it a Boolean dynamics, Fig. 4 (left) giving an example of attractor with final states (in bold red) different from the initial conditions. The characteristics of the associated interaction matrix W are:
- – in 1948, Delbrück [21] conjectured that positive loops in the interaction graph of a regulatory network was a necessary condition for cell differentiation, i.e. for the existence of multiple attractors of the genes expression; this conjecture has been written in a good mathematical context by Thomas in 1980 [22]; we have proved above that the positive loops were related to the observation of multiple attractors, which definitively gives to the positive loops another signification than to the negative ones, more related to the stability of the system (like in the classical Watt regulator, well known in cybernetics, cf. Fig. 5);
- – in 1992, Kauffman [20] conjectured that the mean number of attractors for a Boolean genetic network with n genes and with connectivity 2 was of the order of √n (see Theorem 4 above). This conjecture is now supported by real observations: we have about 35 000 genes in the human genome and about 200 different tissues, which can be considered as different attractors of the same dynamics. For Arabidopsis thaliana, there is
different tissues [10] and for the Cro operon of the phage , , (lytic and lysogenic attractors) [23,24], with Boolean [25] or discrete multi-level [26] models.

Interaction graph of the flowering regulatory network of Arabidopsis thaliana (right) and an attractor of its Boolean dynamics (left).

The Watt regulator, the prototype of negative regulatory loop in cybernetics.
3.2 The gastrulation regulatory network
If we consider the regulatory network ruling the gastrulation in Drosophila (cf. Fig. 6 and [28]), it is easy to check that

The gastrulation regulatory network (after [27]) (ADK and NDK are respectively the adenylate kinase and the nucleotide diphosphate kinase).

The Waddington chreode or morphogenetic landscape.
3.3 The phage μ lytic-lysogenic attractor [32]
If we consider the operon governing the expression of the phage μ, we obtain the graph given in Fig. 8. It is interesting to notice that

The phage μ operon.
3.4 The mnesic ‘opernet’
We will call in the following mnesic ‘opernet’ the system obtained by merging the genetic (operon) and epigenetic parts (metabolic net) of the system ruling the cotyledonary buds growth (cf. [33,34] and Fig. 9).

The mnesic opernet.
The genetic part of the mnesic opernet has been modelled by a Boolean system [34]:
- – variable PA (resp PB) represents pricking treatment (or any other stress action) on the side A (resp B) of the plant and its value is 1 if treatment has been done, and 0 if not;
- – variable SA (resp SB) represents the discrete part of the operon; we suppose that it contributes to mobilize a morphogenetic cotyledonary material RA (resp RB) responsible for the growth of the apex and of the cotyledonary bud A (resp B).
We will suppose in the following that the variable RA (resp RB) representing the concentration of R on side A (resp B) is continuous and that its velocity dRA/dt (resp dRB/dt) is ruled by a differential system containing a three-switch between the continuous variables T (apex growth metabolites concentration), A and B (cotyledonary buds A and B growth metabolites concentrations on respectively A and B side) (see Fig. 10).

Interaction graph of the epigenetic part as a continuous differential system.
Graph G of Fig. 9 is such that
Then we can write the differential system governing the continuous variables T,A,B,RA and RB as follows:
The two first equations correspond to the dynamics of RA (resp RB) whose concentration derivative at time tdRA(t)/dt (resp dRB(t)/dt) results from an auto-catalytic term σRA(t) (resp σRB(t)) diminished by the term (−kT(t)−4kA(t)/5−kB(t)/5) RA(t) (resp (−kT(t)−4kB(t)/5−kA(t)/5) RB(t)) expressing the inhibition by T, A and B, plus a production of growth metabolites term denoted by F(RA(t)) (T(t)+4A(t)/5+B(t))/5)/(T(t)+A(t)+B(t)) (resp F(RB(t))(T(t)+A(t)/5+4B(t)/5)/(T(t)+A(t)+B(t))), by supposing that the RA (resp RB) consumption is competitively inhibited by the bud growth on its side A (resp B) and by the bud growth on the other side and by considering that Kms and
The equations for the apex and cotyledonary buds growth metabolites concentrations just express that their production comes from RA and RB, with a competitive inhibition by the other sources of growth, plus a linear degradation term.
We suppose now for interpreting a minima the experimental results given in Tables 1 and 2 of [33] that F is an allosteric function of order 4 having two successive inflection points (involving that the protein catalysing the production of apex and buds growth metabolites has four catalytic subunits) as allowed by the Monod–Wyman–Changeux equation (see [35] and Fig. 11).

Allosteric function F, which presents two successive inflection points.
If the value of F′ verifies F(r)−rF′(0)<0 and ν<F(r)/r−F′(r)<2 (krF′(r))1/2−ν, then the differential system above possesses at most 16 stationary states, whose 0 is a stable focus, two (respectively (r,0,σr/(ν+kr)) and (0,r,σr/(ν+kr)) are unstable focuses surrounded by limit cycles α and β, and C(r,r,2σr/(ν+kr)) is an attractor (either stable focus, or limit cycle) as shown in Fig. 12, by supposing known the dynamics of the inhibitory three-switch [36] between T, A and B.

Experimental perturbations in the state space (RA,RB,T).
More generally, if we have an inhibitory n-switch between the Ai's verifying:
The possible configurations of experimental perturbations as reported in Table 1 below and in Tables 1 and 2 of [33] give different trajectories after perturbations, as explained in the following. If nA (resp nB) represents the number of seedlings beginning to elongate on the side of bud A (resp B), then g=(nA−nB)/n, where n=nA+nB, is an asymmetry growth index.
Experimental results about domination growth after decapitation
Pricking treatment | g | Domination |
(1) Non-pricked control | 0.02 | A = B |
(2) 4A with decapitation at onset daylight (dod) | 0.35±0.02 | A<B |
(3) 4A with decapitation at midday (dm) | 0.08±0.15 | A=B |
(4) 1A (dod) | 0.01 | A=B |
(5) 4B (dod) | −0.35 | A>B |
(6) 1A (1 h) 4B (dod) | 0.39 | A<B |
(7) 2A (dod) | 0.06 | A=B |
(8) 2A (1 h) 2A/2B (dod) | 0.32 | A<B |
(9) 2A (1 h) 2A/2B (3 h) 2A/2B (dod) | 0.05 | A=B |
(10) 2A (1 h) 2A/2B (3 h) 2A/2B (5 h) 2A/2B (dod) | 0.34 | A<B |
We can make in Fig. 12 above the following observations.
- (1) If the initial conditions are inside the basin of stability of the attractor C, near C, then the whole trajectory without perturbation lies in this basin and tends to the attractor C (where RA≈RB) as time t is tending to infinity. After decapitation without primitive pricking (non-pricked control), A and B buds have the same chance to begin to grow up, then to inhibit the bud growth on the opposite side, hence g≈0.
- (2–3) If we do four pricking treatments on A, then from initial conditions near C we observe a shift to the basin of the limit cycle β and following the decapitation time, we get a B domination if decapitation is done at the onset daylight (dod) and a symmetry in buds growth if it occurs at midday (dm) (see Fig. 12): it is due to the fact that RA and RB after decapitation lies in the basin of the limit cycle β where RB>RA in the first case and are in the basin of the stable focus 0 in the second case.
- (4) If the system is starting near C and if one pricking 1A is done on cotyledon A (resp B), then we suppose that the perturbation results in a decrease of intensity w of RA (resp RB) and in a decrease of intensity w/2 of RB (resp RA) at time tP of pricking; then the variables RA and RB remain in the basin of C.
- (5) If we are doing four pricking treatments on side B, then the point representing the system leaves the basin of C with a decrease of 4w in RA and 2w in RB to go to the basin of the limit cycle β.
- (6) If we do one pricking on side A and 1 h after four prickings on side B, then the system goes in the basin of the limit cycle α and for opposite reasons than for (5) above, RB<RA and A dominates B after decapitation.
- (7) If we do two prickings on side A, then the system remains in the basin of C.
- (8) If we do two prickings on side A and 1 h after two prickings on each side, then the system goes in the basin of the limit cycle β.
- (9) If we do two prickings on side A, 1 h after two prickings on each side and 3 h after two prickings on each side, then the system returns in the basin of C due to the special shape of the trajectory (8) going to the limit cycle β passing near the frontiers of the basin of C.
- (10) If we do two prickings on side A, 1 h after two prickings on each side, 3 h after two prickings on each side and 5 h after anew two prickings on each side, the system goes in the basin of C like in (9) after 4 h, then 5 h after it goes from C to the basin of the limit cycle β.
4 Conclusion
An important first conclusion we have to make explicit in this paper concerns the relationship between the number F of fixed points and the number S of interaction circuits of the interaction matrix W: the problem is in fact to find the best upper bound for F for a given interaction matrix W. This question is related to the famous 16th Hilbert's problem, whose one of the aim is to give an efficient upper bound for the number of limit cycles of a polynomial differential system. Let us summarize the role of the architecture of positive and negative circuits of W on the occurrence of multiple stationary behaviours, as obtained above: if the number of nodes and the number of arcs are the same, there is only one isolated interaction circuit (S=1) in W and either this circuit is negative and the lowest bound (0) for F is reached, or this circuit is positive and the upper bound (21) for F is reached. If the number of nodes is n and the number of arcs is n+1, there is two interaction circuits (S=2) with the following structure: if both circuits are negative, F=0; if there is a positive circuit and a negative circuit disjoint, F=0; if there is a positive circuit intersecting a negative circuit, F=1; if there is a positive circuit intersecting a positive circuit, F=1; if there is two disjoint positive circuits, F=22. If, more generally, the number S of interaction circuits of W is m, then: if all circuits are negative, F=0; if all circuits are positive, 2⩽F⩽2m and if all circuits are positive and disjoint, F=2m. An interesting open problem is now to make exhaustive the determination of F and S and in particular to find the circumstances (related to the circuits structure) in which we can relate the number of intersecting and isolated circuits to F. A conjecture we could make is that F=2c, where c is the number of not-singleton connected components of the interaction graph G having at least one positive loop: it holds for the lactose, Arabidopsis, phage μ and gastrulation regulatory systems. The approach for solving this open problem could consist in finding coherent relationships between analogous properties discovered independently for continuous and Boolean versions of regulatory networks [4–9].
The second conclusion concerns the practical use of the presented results; a geneticist can for example exploit the minimality results in the following sense: we have shown in the paper that it would be possible to characterize the minimal interaction matrices having certain state vectors as fixed points. The determination of these matrices is not unique, but permits to focus on a certain important equivalence class to which the expected matrix has to belong. This considerably restricts the choice of the possible interaction matrices compatible with observed fixed points, when it is impossible to directly get from experiments all interaction coefficients, but when it is only possible to observe the phenomenology of fixed points or limit cycles. This corresponds in genetics to the observation of stationary expression behaviours (for example from bio-array imaging) without experimental measure of the inhibitory and activatory coefficients of repressors and promoters. The possibility to obtain (even in an equivalence class) a sketch of the interaction matrix permits to construct (by randomising in a Bayesian way the unknown coefficients of
Acknowledgements
We have done this work thanks to the support of the National Network for Technology Research RNTS ‘Technologies for Health’ from the French Ministry of Research.