## 1 Introduction

In living organisms, many proteins are transcription factors: they can bind to DNA and regulate the transcription of specific genes, i.e. the synthesis of RNA from coding regions of chromosomal DNA. This regulation of transcription is a very complex mechanism, which can involve up to dozens of genes, and other factors as well. Furthermore, regulation occurs during the full process of gene expression. In all cases, one will say that a gene A activates (resp. inhibits) a gene B when A produces a protein that has a positive (resp. negative) effect on the expression of gene B. If several genes are involved in a given biological system, they form a gene network from which one can draw an interaction graph G. In mathematical terms, G is a finite oriented graph, the edges of which are endowed with a sign: the vertices are the genes, and a positive (resp. negative) edge from j to i means that j activates (resp. inhibits) i. Note that a gene can activate or inhibit itself, i.e. G can have edges that end where they start. Furthermore, depending on the concentrations of the proteins in the system, the effect of j on i can be positive, negative, or absent. In other words, G is a function of the concentrations.

In general, very little is known about the strength of the interactions between genes. One is thus faced with the following difficult problem: which dynamical properties of a gene network can be inferred from the topology of its interaction graph (despite the lack of quantitative information)? Several methods have been used to tackle this question. One of them is numerical simulation, which requires choosing kinetic parameters in a realistic way. Another method is to study the statistical properties of gene networks, by comparing their interaction graphs with random ones. One can also try to decompose a given graph into submodules of biological significance. Finally, some authors have focused their attention on special motifs, i.e. subgraphs with simple topology, for instance, those involving few vertices that are overrepresented in gene networks [1].

In this paper we shall study circuits. A circuit C in an interaction graph G is a sequence ${e}_{1},\dots ,{e}_{k}$ of edges such that the end point of ${e}_{i}$, $i=1,\dots ,k-1$ (resp. ${e}_{k}$) is the origin of ${e}_{i+1}$ (resp. ${e}_{1}$), each vertex of G occurring at most once in C. The sign of a circuit is the product of the signs of its edges. When gene regulation was discovered, it was soon noticed that circuits (or ‘feedback loops’) are often present in gene networks (or at least in mixed networks, with interactions between genes, proteins and metabolites), and that their biological role depends on their sign. For instance, if $G=C$ consists of a single positive circuit, the network can have two possible stable stationary states. For instance, when G looks as follows:

However, when the interaction graph G contains a positive circuit C, there is no reason to expect that C will govern the dynamic behaviour of the underlying network; this depends on the strength of the interactions. Thomas [3] had the idea to turn the statement another way: positive circuits are necessary if not sufficient for multistationarity. He proposed the following:

Thomas' rule: Assume a gene network has several non-degenerate stationary states. Then its interaction graph contains, somewhere in phase space, a positive circuit.

Here, the expression “somewhere in phase space” refers to the fact, mentioned earlier, that the concentrations of proteins have to be given for the graph to be defined.

This rule of Thomas, if true, can be quite useful for geneticists. For instance, if one knows that a given biological system can differentiate, one will be led to search for a positive circuit relating the genes involved. It is therefore worthwhile to decide how general this rule is. One way to check it is the following. First, describe a gene network in a formal mathematical way (a model). Then, within this formalism, make sense of all the terms figuring in Thomas' rule (non-degenerate stationary state, interaction graph, phase space...). Thomas' rule then becomes, in this context, a precise mathematical statement (a conjecture), which one can try to prove (or to contradict) logically.

For instance, Thomas and Kaufman phrased Thomas' rule [3] as a precise mathematical conjecture by using a differential model of gene networks [4]. This conjecture was proved under additional assumptions in [5–7], and in general in [8]. Later, Thomas' rule was checked for Boolean models [9]. Many more models of gene networks have been proposed (see [10] for a thorough survey). In this paper, we shall show that Thomas' rule is true for five different types of models of gene networks: the Boolean, differential, differential with decay, piecewise-linear and multivalued discrete models. Of course, knowing that the rule is true for one type of models does not imply automatically that it is true for another one. Still, our arguments will exhibit interesting connections between several ways of modelling. When the spontaneous decay of all proteins is taken into account, Thomas' rule happens to be more robust, and this allows us to study the piecewise-linear case by approximation. And, in some cases, the piecewise-linear models can in turn be described in discrete terms. We conclude the paper by discussing open problems. Among them is whether Thomas' rule remains valid for other biological networks, or when the stochastic nature of gene regulation is taken into account.

## 2 Boolean models

Let $n\u2a7e1$ be an integer and $\Omega ={\{0,1\}}^{n}$, the set of strings of n letters in the alphabet $\{0,1\}$. Consider a map

$$F=({F}_{i})\phantom{\rule{0.2em}{0ex}}:\Omega \to \Omega $$ |

To every x in Ω we attach an interaction graph $G(x)$ which is described as follows. Fix $j\in \{1,\dots ,n\}$ and let $y\in \Omega $ be defined by:

$${y}_{j}=1-{x}_{j}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{y}_{k}={x}_{k}\phantom{\rule{1em}{0ex}}\text{if}k\ne j$$ |

A stationary state of the network is a fixed point of F, i.e. a point $x\in \Omega $ such that $F(x)=x$. Part (2) of the following theorem says that Thomas' rule is true for Boolean models.

#### Theorem 1

## 3 Differential models

Let $n\u2a7e1$ be an integer and $\Omega ={\mathbf{R}}^{n}$ the standard real vector space of dimension n. Consider a differential map:

$$F=({F}_{i})\phantom{\rule{0.2em}{0ex}}:\Omega \to \Omega $$ |

$$\frac{\mathrm{d}x}{\mathrm{d}t}=F(x)$$ | (1) |

Given $x\in \Omega $, we define an interaction graph $G(x)$ as follows. Its set of vertices is $\{1,\dots ,n\}$ and there is a positive (resp. negative) edge from j to i when the partial derivative $\frac{\partial {F}_{i}}{\partial {x}_{j}}(x)$ is positive (resp. negative). [To illustrate this, assume $\frac{\partial {F}_{i}}{\partial {x}_{j}}(x)>0$. If we increase the concentration of protein j in state x, the number ${F}_{i}(x)$ will increase, and the production of i will accelerate. In other words, j is an activation of i in state x.]

From (1), we see that a stationary state of the network is a zero of F, i.e. x ∈Ω such that $F(x)=0$. This zero is called non-degenerate when the determinant $\mathrm{det}(\frac{\partial {F}_{i}}{\partial {x}_{j}}(x))$ of the Jacobian matrix at x is different from zero. Thomas' rule is true for differentiable models.

#### Theorem 2 [8]

Assume that F has at least two non-degenerate zeroes. Then there exists$x\in \Omega $such that$G(x)$contains a positive circuit.

For previous results, see [5–7].

## 4 Differentiable models with decay

Since concentrations cannot be negative, one would like to get a version of Theorem 2 where $F(x)$ is only defined for those $x=({x}_{i})$ such that ${x}_{i}\u2a7e0$, $i=1,\dots ,n$. But it turns out that it is not true as stated with this restriction [8 (§3.5)]. However, a more realistic modelling of gene networks consists in taking into account that the concentration of every protein is submitted to a spontaneous decay, due to degradation and to the growth of cells. We are thus led to the following model.

For every $i=1,\dots ,n$ let ${\Omega}_{i}\subset \mathbf{R}$ be a real interval (i.e. ${\Omega}_{i}=[{a}_{i},{b}_{i}],]{a}_{i},{b}_{i}],[{a}_{i},{b}_{i}[$ or $]{a}_{i},{b}_{i}[$, with $-\infty \u2a7d{a}_{i}$ and ${b}_{i}\u2a7d+\infty $). On the product $\Omega ={\prod}_{i}{\Omega}_{i}\subset {\mathbf{R}}^{n}$, consider a differentiable map

$$F=({F}_{i})\phantom{\rule{0.2em}{0ex}}:\Omega \to {\mathbf{R}}^{n}$$ |

$$\frac{\mathrm{d}{x}_{i}}{\mathrm{d}t}={F}_{i}(x)-{\gamma}_{i}{x}_{i},\phantom{\rule{1em}{0ex}}i=1,\dots ,n$$ | (2) |

For every $x\in \Omega $, let $G(x)$ be the interaction graph defined from the signs of the partial derivatives $\frac{\partial {F}_{i}}{\partial {x}_{j}}(x)$ as in §2 above.

#### Theorem 3

Assume that there exists two points$x\ne y$in Ω such that

$$|{F}_{i}(x)-{\gamma}_{i}{x}_{i}-{F}_{i}(y)+{\gamma}_{i}{y}_{i}|<{\gamma}_{i}|{x}_{i}-{y}_{i}|$$ |

#### Remarks

- (1)
A stationary state is a point $x\in \Omega $ such that, for all $i=1,\dots ,n$,
$${F}_{i}(x)-{\gamma}_{i}{x}_{i}=0$$ If $x\ne y$ are two stationary states of (x) we have, for all $i=1,\dots ,n$

therefore the hypotheses of Theorem 3 are satisfied. In other words, Thomas' rule is again true, for arbitrary Ω, when we take the spontaneous decay of proteins into account. It is also more ‘robust’, since it remains valid for two states $x\ne y$ that are almost stationary.$${F}_{i}(x)-{\gamma}_{i}{x}_{i}={F}_{i}(y)-{\gamma}_{i}{y}_{i}=0$$ - (2) In the conclusion of Theorem 3, we can assume that ${z}_{i}$ = ${x}_{i}$ whenever ${x}_{i}={y}_{i}$. To see that, let I be the set of indices i such that ${x}_{i}={y}_{i}$, and apply Theorem 3 to the restriction of the functions ${F}_{i}$, $i\notin I$, to the linear subspace of those z $\in {\mathrm{R}}^{n}$ such that ${z}_{i}={x}_{i}$ when $i\in I$. A similar remark can be made in Theorems 1, 2 and 4.

## 5 Piecewise-linear models

Let $\Omega ={\prod}_{i}{\Omega}_{i}$ be as in §4. Another model for gene networks [13–15] is given by the system of equations

$$\frac{\mathrm{d}{x}_{i}}{\mathrm{d}t}={F}_{i}(x)-{\gamma}_{i}{x}_{i}$$ |

More precisely, for every real number θ, define the step function $s(x,\theta )$ from **R** to subsets of **R** by:

$$s(x,\theta )=\{\begin{array}{cc}1\hfill & \mathrm{if}\phantom{\rule{0.25em}{0ex}}x>\theta \hfill \\ ]0,1[\hfill & \mathrm{if}\phantom{\rule{0.25em}{0ex}}x=\theta \hfill \\ 0\hfill & \mathrm{if}\phantom{\rule{0.25em}{0ex}}x<\theta \hfill \end{array}$$ |

For every $j=1,\dots ,n$, choose finitely many distinct real thresholds ${\theta}_{j}^{k}$ in the interior of ${\Omega}_{j}$, $k=1,\dots ,{m}_{j}$. For each $i=1,\dots ,n$, fix a real polynomial ${P}_{i}({T}_{j}^{k})$ in ${\sum}_{j=1}^{n}{m}_{j}$ variables and, for every x in Ω, let:

$${F}_{i}(x)={P}_{i}\left(s({x}_{j},{\theta}_{j}^{k})\right)$$ | (3) |

By definition, ${F}_{i}(x)$ is a subset of **R**, reduced to a single point if ${x}_{j}\ne {\theta}_{j}^{k}$ for all j and k.

Now we consider the system of piecewise linear differential inclusions:

$$\frac{\mathrm{d}{x}_{i}}{\mathrm{d}t}\in {F}_{i}(x)-{\gamma}_{i}{x}_{i},\phantom{\rule{1em}{0ex}}i=1,\dots ,n$$ | (4) |

A stationary state of (4) is a point x in Ω such that 0 lies in (${F}_{i}(x)-{\gamma}_{i}{x}_{i}$). For every x in Ω, we define an interaction graph $G(x)$ as follows. Its set of vertices is $\{1,\dots ,n\}$ and there is a positive (resp. negative) edge from j to i when ${x}_{j}={\theta}_{j}^{k}$ is a threshold and the value of the partial derivative $\partial {P}_{i}/\partial {T}_{j}^{k}$ is positive (resp. negative) at the some point in the set ($s({x}_{a},{\theta}_{a}^{b})$). Note that in the case considered in [16] the interaction graph of §1 of op. cit. is the superposition of all the graphs $G(x),x\in \Omega $.

#### Theorem 4

Assume that(4)has several stationary states. Then there exists$x\in \Omega $such that$G(x)$contains a positive circuit.

## 6 Discrete models

In [16] (Theorem 1) and [17] (Theorem 2), it is shown that the stationary states of some piecewise linear models can be described by fixed points of a map $\Omega \to \Omega $, where Ω is a product of n finite sets ${\Omega}_{1},\dots ,{\Omega}_{n}$, with ${\Omega}_{i}$ of order ${n}_{i}+1$, where ${n}_{i}$ is the number of threshold values of the variable ${x}_{i}$. Theorem 4 gives therefore a proof that Thomas' rule holds for the fixed points of these discrete models.

For example, assume that ${\Omega}_{i}=\{0,1\}$ for all $i=1,\dots ,n$ and let $F=({F}_{i})\phantom{\rule{0.2em}{0ex}}:\Omega \to \Omega $ be any map as in Theorem 1 above. Consider the system of piecewise-linear differential equations:

$$\frac{\mathrm{d}{x}_{i}}{\mathrm{d}t}={\tilde{F}}_{i}(x)-{x}_{i},\phantom{\rule{1em}{0ex}}i=1,\dots ,n$$ | (5) |

$$d\phantom{\rule{0.2em}{0ex}}:\tilde{\Omega}\to \Omega $$ |

$$d{(x)}_{i}=(1+\mathrm{sign}({x}_{i}))/2$$ |

$${\tilde{F}}_{i}(x)=\{\begin{array}{cc}1\hfill & \mathrm{if}\phantom{\rule{0.25em}{0ex}}{F}_{i}(d(x))=1\hfill \\ -1\hfill & \mathrm{if}\phantom{\rule{0.25em}{0ex}}{F}_{i}(d(x))=0\hfill \end{array}$$ |

$$d(\tilde{F}(x))=F(d(x))$$ |

Then one can check that, for every $x\in \tilde{\Omega}$, if y lies in the closure in ${\mathbf{R}}^{n}$ of the component of $\tilde{\Omega}$ containing x, the interaction graph $G(y)$ defined from $\tilde{F}$ as in §4 is contained in the interaction graph $G(d(x))$ of §1. Furthermore, according to Theorem 1 [16], if $d(x)$ is a fixed point of F, the point $x\in \tilde{\Omega}$ is an ‘asymptotically stable’ steady state of (5). From Theorem 4, we thus get a new (and quite indirect!) proof of Theorem 1(2) under the above hypothesis. It would remain to decide which discrete models, and in particular which Boolean models, can be obtained from the piecewise-linear models considered in [16,17].

## 7 Concluding remarks

The results above give support to the validity of Thomas' rule. They do not ‘prove’ any kind of ‘biological law’, since they depend on the quality of the models we used for describing gene networks, and these are, clearly, gross simplifications of the biological reality. For instance, we did not consider the role of chromatin conformation in the regulation of transcription. Neither did we include any discussion of the alternate splicing phenomenon. Therefore, it might be worth checking this rule in new and more refined setups. Let us mention a few possible extensions of the results presented here.

Gene networks do not appear in isolation. They are usually coupled with other biological networks, like the metabolic networks, and those involving interactions between proteins. People have tried to describe these mixed networks in a single picture (see for example [18] or [19]), but this is not so easy, since edges in these enlarged graphs do not have the same meaning as in the case of gene regulation. On the other hand, one should notice that an enlarged mixed network may well contain a positive circuit that is not visible in any smaller pure one, as soon as its vertices are of different natures.

Thomas' rule is of course valid for any system that can be modelled by one of the five methods described above. However, one has to be careful that the interaction graphs defined in these models by means of a mathematical recipe need not have an obvious intuitive meaning. In particular, they might not be simply related to the usual way of representing these systems (as noted in [20], Remark 5)), and a positive circuit needs not be visible in the traditional graphic representation. For instance, let us consider a chemical system containing, among others, two compounds X and Y, and a reversible reaction:

$$\mathrm{d}[X]/\mathrm{d}t=a[Y]+\text{other terms}$$ |

$$\mathrm{d}[Y]/\mathrm{d}t=b[X]+\text{other terms}$$ |

Another way to pursue our discussion would be to take into account the location of the proteins. Since the famous work of Turing in the 1950s, many models of development have considered diffusion phenomena. It might be worth noting, though, that the scenario of spatial differentiation proposed by Turing and its followers requires that at least one of the chemicals be a self-activator. In that sense, there is still a positive circuit involved, and if we delete the diffusion terms in Turing's equations we get back to the situation described in Theorem 1. I do not know if diffusion in space can lead to differentiation in the absence of any positive circuit in the appropriate interaction graph.

A third direction is the following. All the models presented above are deterministic. Now, many recent works on gene expression insist on the importance of stochastic effects. These are manifest in the variation of expression levels from one cell to another. This forces us to view the binding of a protein to DNA, and the whole gene regulation, as a stochastic event. In [23], Gillespie proposed to represent the elementary chemical reactions in a gene network as a discrete jump Markov process, and, when there are enough copies of each protein, to approximate the chemical master equation by Fokker–Planck stochastic differential equations. It would be very interesting to decide if a variant of Thomas' rule still remains true in such a context. It is known that introducing stochasticity in a deterministic model allows for occasional switching between different stationary states (and such switches have been observed experimentally, see [24] for a survey). But can Gillespie's model lead to completely new scenarios of differentiation (violating Thomas' rule)?

Finally, one can also seek upper bounds for the number of possible stationary states in a gene network [11]. Such a bound is obtained in [25] when gene networks are modelled by means of IN and OR networks. Similar upper bounds cannot be valid for an equation like (1), since, obviously, the number of stationary states depends on the complexity of the function F (e.g., its degree when F is algebraic). For a general result on the complexity of gene networks, we refer to [26].

## Acknowledgements

The author wishes to thank H. De Jong, J. Demongeot, M. Kaufman, N. Le Novère, O. Radulescu, E. Rémy, P. Ruet, D. Thieffry, R. Thomas, and A. Wagner for helpful discussions and comments.

## Appendix A Proof of Theorem 3

We proceed by contradiction and assume that none of the graphs $G(z),z\in \Omega $, contains a positive circuit. Then, according to [8] (Lemma 2i), all the principal minors of the Jacobian matrix of $(-{F}_{i})$ at every point z are nonnegative. Fix a constant $\varepsilon >0$. By [8] (4), we conclude that all the principal minors of the Jacobian matrix of ($-{F}_{i}(x)+\varepsilon {x}_{i})$ are positive in Ω. By the univalence theorem of Gale–Nikaido ([27,28 (p. 20a)]), this implies that the map $(-{F}_{i}(x)+\varepsilon {x}_{i})$ is a P-function on Ω, i.e., when $x\ne y$ are two points in Ω, there exists $k\in \{1,\dots ,n\}$ such that:

$$({x}_{k}-{y}_{k})(-{F}_{k}(x)+\varepsilon {x}_{k}+{F}_{k}(y)-\varepsilon {y}_{k})>0$$ |

$$({x}_{k}-{y}_{k})(-{F}_{k}(x)+{F}_{k}(y))\u2a7e0$$ |

$$({x}_{k}-{y}_{k})({F}_{k}(x)-{\gamma}_{k}{x}_{k}-{F}_{k}(y)+{\gamma}_{k}{y}_{k})\u2a7e{\gamma}_{k}{({x}_{k}-{y}_{k})}^{2}$$ |

$$|{F}_{k}(x)-{\gamma}_{k}{x}_{k}-{F}_{k}(y)+{\gamma}_{k}{y}_{k}|\u2a7e{\gamma}_{k}|{x}_{k}-{y}_{k}|$$ |

## Appendix B Proof of Theorem 4

Following a suggestion of J.-L. Giavitto, we use an argument of approximation. For every θ∈ **R** and every integer m ⩾0, choose a differentiable function ${s}_{m}(x,\theta )$ on **R** such that ${s}_{m}(x,\theta )=s(x,\theta )$ when $|x-\theta |\u2a7e\frac{1}{m}$, ${s}_{m}(x,\theta )$ is strictly increasing when |x $-\theta |<\frac{1}{m}$, and, given x ∈ **R** and u $\in s(x,\theta )$, for every $\varepsilon >0$, there exists ${m}_{0}$ such that, if $m\u2a7e{m}_{0}$, there is $\xi \in \mathbf{R}$ with:

$$|\xi -x|<\varepsilon $$ |

$$|{s}_{m}(\xi ,\theta )-u|<\varepsilon $$ |

$$|{\xi}_{i}-{x}_{i}|<\varepsilon $$ |

$$|{F}_{i}(\xi ,m)-{u}_{i}|<\varepsilon $$ |

$${F}_{i}(\xi ,m)={P}_{i}\left({s}_{m}({x}_{j},{\theta}_{j}^{k})\right)$$ |

Assume now that the system (4) has two non-degenerate stationary states $x\ne y$ in Ω. The assertion $0\in ({F}_{i}(x)-{\gamma}_{i}{x}_{i})$ means that ${u}_{i}\in {F}_{i}(x)$, where ${u}_{i}={\gamma}_{i}{x}_{i}$. Similarly, ${v}_{i}\in {F}_{i}(y)$ with ${v}_{i}={\gamma}_{i}{y}_{i}$. For every $\varepsilon >0$, and m big enough, choose ξ and η in Ω such that, for all $i=1,\dots ,n$:

$$|{\xi}_{i}-{x}_{i}|<\varepsilon ,\phantom{\rule{2em}{0ex}}|{\eta}_{i}-{y}_{i}|<\varepsilon ,$$ |

$$|{F}_{i}(\xi ,m)-{u}_{i}|<\varepsilon ,\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}|{F}_{i}(\eta ,m)-{v}_{i}|<\varepsilon $$ |

From these inequalities we conclude that:

$$|{F}_{i}(\xi ,m)-{\gamma}_{i}{\xi}_{i}-{F}_{i}(\eta ,m)+{\gamma}_{i}{\eta}_{i}|<|{u}_{i}-{\gamma}_{i}{\xi}_{i}-{v}_{i}+{\gamma}_{i}{\eta}_{i}|+2\varepsilon <2{\gamma}_{i}\varepsilon +2\varepsilon $$ |

$$2{\gamma}_{i}\varepsilon +2\varepsilon <{\gamma}_{i}|{\xi}_{i}-{\eta}_{i}|$$ |

It follows that $({F}_{i}(\cdot ,m))$ satisfies the assumption of Theorem 3, hence its interaction graph contains a positive circuit somewhere in phase space. It remains to notice that, for every z ∈Ω, when m is big enough, the interaction graph of (${F}_{i}(.,m)$) at z is the same as the one of F at some point in Ω.