Theoretical Study of the 1 Self-Regulating Gene in the Modiﬁed Wagner Model

: Predicting how a genetic change affects a given character is a major challenge in biology, and being able to tackle this problem relies on our ability to develop realistic models of gene networks. However, such models are rarely tractable mathematically. In this paper, we propose a mathematical analysis of the sigmoid variant of the Wagner gene-network model. By considering the simplest case, that is, one unique self-regulating gene, we show that numerical simulations are not the only tool available to study such models: theoretical studies can be done too, by mathematical analysis of discrete dynamical systems. It is ﬁrst shown that the particular sigmoid function can be theoretically investigated. Secondly, we provide an illustration on how to apply such investigations in the case of the dynamical system representing the one self-regulating gene.


Introduction
Predicting the effect of a genetic change (differences in the DNA molecule) on a character of interest (which can be related to, e.g., human health, plant and animal production, or evolutionary differences between species) remains a major challenge in biology [1].This is mainly due to the fact that cell physiology is heavily regulated by complex gene networks, which are able to compensate various genetic defects or environmental disturbances.Understanding how DNA variation translates to observable (phenotypic) variation, due for instance to a modification of the conformation of proteins [2] or DNA [3] a is of obvious importance, although very complex.Being able to tackle this complexity relies, among other things, on our ability to develop realistic models of such gene networks.
There exist several kinds of biological interactions that can be modeled as gene networks (e.g., signal transduction, metabolism, or transcription regulation).The literature proposes several modeling frameworks for each of them based on various biological hypotheses and different time scales [4][5][6].An interesting subset of gene network models does not aim at predicting the behavior of a specific group of identified genes in an organism, but are rather used as a general abstraction of a gene network, in order to study their evolutionary properties in individual-based simulations [7,8].Although naive in terms of biological hypotheses, these models are particularly important because they could contribute to unifying systems biology and evolutionary genetics.
The most popular framework in this field is the so-called "Wagner Model" (after [7,9]).This model is an abstraction of the interactions between transcription factors (genes that regulate the expression of other genes).The structure of the n-gene network is encoded into a n × n matrix (W in the original model, thereafter noted M), which is constant for an organism, and the status of the gene network (the expression of the n genes of the network) is encoded into a vector S of size n.The dynamics of the gene network over a finite number of discrete time steps is modeled as a system of n linear difference equations, often written as S t+1 = F(MS t ), where F(x) = [ f (x 1 ), f (x 2 ), ..., f (x n )], in which f is a scaling function.The n 2 elements of matrix M are real numbers representing the influence of gene j on gene i (M ij < 0 for repression, M ij > 0 for up-regulation, and M ij = 0 for the absence of interaction).Self regulation is possible (and realistic), i.e., S ii can be different from 0. The purpose of the scaling function f is to ensure that gene expressions S remain in their domain of definition.The gene network starts from an initial vector of gene expressions S 0 , generally set at an arbitrary level.
In the original setting from [9], gene expressions were discrete and could take only two values, -1 for no expression and 1 for an expressed gene, f (x) was thus a step function ( f (x < 0) = −1; f (0) = 0; f (x > 0) = 1).Alternative parameterization include e.g.scaling between 0 and 1 [10].An overview of the diversity of similar settings is provided in [11].Recent implementations of the model consider continuous gene expressions, and the step function was turned into a sigmoid, scaling gene expressions between -1 and 1 ( f (x) = 2/(1 + e −x ) − 1, [12], see [13] for a mathematical analysis).For more realism, the sigmoid function can also be further modified to ensure that genes are only weakly expressed in absence of regulators, by considering that f (0) = a < 1/2, as in [14], and which is the model studied below.
The main purpose of such models is to ensure that any combination M, S 0 can be solved computationally (e.g., as the state of the network S T after T timesteps) within a predictable (most of the time, constant) amount of time.This is of major importance in individual-based computer simulations or other numerical studies in which the network structure M can mutate and evolve over time.The lack of mathematical tractability remains, however, problematic, as it makes it difficult to compare simulation results with classical predictions from population and quantitative genetics models [15,16].This is why, in this article, a mathematical analysis of the sigmoid variant of the Wagner gene-network model is presented, focusing on the simplest case, that is, one unique self-regulating gene.It is demonstrated that, in addition to classical numerical simulations, this model can be theoretically studied too, by the means of mathematical analysis of discrete dynamical systems.It is first shown that the particular sigmoid function, usually studied within this model, can be theoretically investigated, and such investigations are secondly partially applied to the dynamical system representing the one self-regulating gene.Finally, some ways to extend the analysis to the multiple gene case are sketched.
The remainder of this article is structured as follows.In the next section, the parameterized sigmoid function is deeply studied by the means of mathematical analysis.Effects of parameter changes on its shape are discussed too, and a systematic investigation of its fixed points is finally provided.Section 3, for its part, focuses on the discrete dynamical system of the modified Wagner model.The one self-regulating gene model is next partially studied, in a particular situation and then in its most general formulation.This article ends by a conclusion section, in which the contribution is summarized and intended future work is outlined.

Studying the sigmoid function
We first provide some rationales about the sigmoid function usually considered in the variant of the Wagner gene-network model studied here [14].They will be used in the next section, when studying the one self-regulating gene system.

Introducing the considered sigmoid
Let us consider a ∈]0, 1[ and the particular sigmoid function defined by: This function, and its parameters, have been chosen in order to have: • a continuous increasing function, • the limit of f a as x approaches negative infinity is 0, • the limit of f a as x approaches infinity is 1, which corresponds to the curve depicted in Figure 1.

Preprints
( We can thus verify that f a is strictly increasing.Additionally, we can reformulate this derivative, to relate it to the logistic map: (3)

About λ and µ parameters
Let us now investigate the two parameters inside f a that both depend on a.
, for its part, has a derivative equal to: whose variation table is as follows:

+∞ +∞
Let us thus remark that λµ = One of the most important elements to study, when investigating the Wagner model, is the existence and meaning of fixed points.We will show that the fixed points of the one self-regulating gene model are related to those of the modified sigmoid function f a , which are studied hereafter.
We can first remark that, as ∀a ∈ [0, 1], f a is continuous and such that f a ([0, 1]) ⊂ f a ([0, 1]), it is thus a function mapping a compact convex set to itself.Due to Brouwer's fixed-point theorem, f a has at least one fixed point in [0, 1].Let us prove that, Proof.Let us first remark that, if f a has a fixed point x, then it satisfies: (5) Moreover, as f a (x) = x on the one hand, and f a (R) =]0, 1[ on the other hand, we necessarily have: Additionally, 0 < x and f a is strictly increasing, so f a (0) < f a (x), which leads to a < x.To sum up, if f a has a fixed point x, then this latter satisfies: Let us consider a fixed point x for f a .So we have: Depending on the sign of λe −2 − 1, H can be either negative on the whole R set, or positive on a bounded interval containing 2 µ .More precisely, , then H 0 on R, and we have lim , then H(x) = 0 has one unique solution, i.e., f a has one unique fixed point.
, then there exist two real numbers x 1 and x 2 such that the following variation table .
Furthermore, H 1 (as H > 0 on this interval).As H (x 1 ) = 0, we can deduce that ) is negative if and only if j(x 1 ) 0.
We now investigate the sign of j.
, it is sufficient to study j on the interval I = 0, 1 2 .
, so j is strictly decreasing on I.The discriminant of the quadratic equation j(x) = 0 being µ(µ − 4) 0, the latter has two solutions µ ± µ(µ − 4) 2µ , which are equal when µ = 4 (i.e. when a = ).Note that only µ − µ(µ − 4) 2µ may belong to I, and that the latter is equal to . We successively have which is positive for each a ∈]0; 1[ and less than 1 2 .Thus µ − µ(µ − 4) 2µ ∈]0; 1/2[ and we have:  Let us show that , we then have a − 1 2 < 0, and so As stated at the beginning of the proof, since x 1 is a root of H, x 1 is a fixed point for f a and thanks to (7).
a < x 1 , and so j(x 1 ) < 0, as j is decreasing with j(a) = 0. To put it in a nutshell, H(x 1 ) > 0.
A numerical simulation based on a dichotomic approach to solve the equation x on [0, 1], leads to the curve depicted in Figure 3.
Let us now establish the following result: Proof.As stated previously (equation ( 3)), . So, we have Thanks to inequation (2) it can be deduced that f a (x) shares then the same sign than 1 − 2 f a on R, and consequently on [0, 1] where the fixed point x is located.Furthermore, To sum up, if a 1 2 , then f a is a contraction mapping such that | f a | is bounded by a.
Let us now consider a < 1 2 , and x = ln(λ) and each term of the right-side product is in [0, 1].Additionally, as 1 − 2 f a (x) 0 ⇐⇒ x ln(λ) µ = x and as x < 1, we can deduce that the sigh of f a , which is the one of 1 − 2 f a (x), is negative for x 1.
From all the material detailed previously, we can thus conclude that, Theorem 2. The unique fixed-point of f a is an attractive one.
Proof.As | f a (x)| < 1, f a is a contraction mapping.By applying the Banach fixed-point theorem, we can deduce again the existence and uniqueness, and in addition the exponential convergence of the dynamical system to this fixed point.
Note that: 1.The unique fixed-point can be found as follows: start with an arbitrary element x 0 in R and define a sequence (x n ) n∈N by x n = f a (x n−1 ), then x n −→ x(a).2. For a 1/2, as | f a (x)| < a, we can deduce that f a is Lipschitz continuous, with a Lipschitz constant equal to a.As a well-known consequence, the convergence of the aforementioned sequence is at least geometric, with a common ratio of a.For a < 1/2 the same conclusion holds but for a constant γ < 1.

The discrete dynamical system under consideration
Let us firstly recall the general model to study gene networks.Let n ∈ N * , M ∈ M n (R) be a square real matrix of size n × n, and let X 0 ∈ R n .We consider the discrete dynamical system: where Note that, most of the times, X 0 = (a, . . ., a) T .We will now focus on its most simplest cases.

A fundamental case: M = (1)
Let us consider first that n = 1 and that the matrix M is the identity: M = (1).The (Σ) system becomes: Let us recall that, if the recurrent sequence u k+1 = f (u k ) converges, then the limit is a fixed point of f .And, from the study of the previous section concerning the sigmoid f a , we know that this fixed-point x(a): 1. exists and is unique, 2. is such that 0 a < x(a) < 1, 3. the convergence speed is geometric, of ratio equal to a (γ).
As f a is increasing, then the sequence (x k ) k∈N is monotonic.Being bounded, as f a outputs are in ]0, 1[, we can conclude that the sequence converges, and thus its limit is x a .Finally, if x 0 < x(a), then the sequence (x n ) n∈N increases to its limit, and otherwise it decreases to its limit x(a), as depicted in Figure 4.

General 1-D case: M = (m)
In that case, (Σ) becomes: Let us introduce g a,m (x) = f a (mx).

Fixed points of g a,m
We first investigate the fixed points of g a,m , for m > 0 : The case m < 0 will be approach later in a different way.As in the fundamental case of Section 2.  x , then λe −2 − 1 is negative and then H 0 over R, and after the computation of the limits of H m as x approaches ±∞, we can deduce the following table of variations (which is independent of m): x , then H m (x) = 0 has one unique solution, i.e., g a,m has one unique fixed point.

•
If a 1 e 2 − 1 , then we a curve similar to the fundamental case for H m .
x As previously, we remark that H (x As a consequence, .
Again as previously, H m , and H m (x 1 ) = 0, so x 1 > 1 µm .Consequently: and so H m (x 1 ) shares the same sign than j m (x 1 ).
Let us study the quadratic polynomial j m (x) on R. Its discriminant ∆(j m ) is equal to µ 2 m 2 − 4µm = µm(µm − 4), and it has the following sign table: can conclude that j m (x 1 ) > 0. So H m (x 1 ) < 0, and H m thus changes twice its sign, respectively in ]−∞, x 1 [ and in x 1 , 2 µm .Thus g m,a has one fixed point in each of these two intervals.The existence of a unique fixed point can be proved without computations, however deeper investigation to establish the nature of this fixed pointed are necessary.
Lemma 1. g a,m (x) has at least one fixed point.
Proof.Suppose for an absurdum that g a,m has no fixed points.Since g a,m is a continuous map and g a,m (0) = a we have Γ(g a,m ) := {(x, y) : g a,m (x) = y} ⊂ {(x, y) : x < y}, that is trivially false since there are x ∈ R such that x > g a,m (x), (for instance x = 2/m).Theorem 3. If m < 0, then g a,m has a unique fixed point.
Proof.Let P ⊂ R 2 be the subset of the fixed points of g a,m , i.e.P = {(x, x) ∈ R 2 : g a,m (x) = x}; since P is a compact set, we can consider his minimum p 1 ∈ P, that we may call the first fixed point of g a,m .In fact P = Γ(g a,m ) ∩ Γ(y = x) ∩ [0, 1] 2 , thus P is compact because is a finite intersection of close subsets contained in a compact set ([0, 1] 2 ).
Suppose for an absurdum to have another fixed point p 2 ; since p 1 is the first fixed point, we have p 1 < p 2 .However g a,m is a decreasing function on [0, 1], in fact g a,m (x) = λµme −µmx (1 + λe −µmx ) 2 , thus we have p 1 = g a,m (p 1 ) > g a,m (p 2 ) = p 2 , a contradiction that proves that p 1 can be the only fixed point of g a,m .
As can be seen, the number of fixed-point of g a,m can be deduced theoretically, according to the values of a and m.At each fixed-point x, we still have to study its attractive property thanks to the value of |g a,m (x)|.Although technical, this latter can be done by using usual methods from the mathematical analysis.

Conclusion and future work
In this article, the objective was to show that gene-network models like the Wagner one, based on the iterations of some discrete dynamical systems defined using a sigmoid function, can be studied too theoretically.Iterations of the sigmoid function has been deeply studied in a first section, emphasizes the fact that the number of fixed points, their approximative location, and their attractivity, can be computed mathematically.Furthermore, we shown that each iteration of this sigmoid function tends to the fixed-point, no matter the initial condition.Elements showing how to extend this study to the complete one self-regulating gene has then be provided in a second section, showing the possibility of such studies.In future work, we intend to extend such investigations to a network of more than one gene.We will first deeply studied the case of 2-4 genes.Such investigations will then be extended to a larger number of genes, introducing qualitative methods that will depend on the shape of the considered matrix.

Figure 2 .
Figure 2. Parameter dependence against a