Type-dependent stochastic Ising model describing the dynamics of a non-symmetric feedback module

We study an alternative approach to model the dynamical behaviors of biological feedback loop, that is, a type-dependent spin system, this class of stochastic models was introduced by Fern\'andez et. al (2009), and are useful since take account to inherent variability of gene expression. We analyze a non-symmetric feedback module being an extension for the repressilator, the first synthetic biological oscillator, invented by Elowitz and Leibler (2000). We consider a mean-field dynamics for a type-dependent Ising model, and then study the empirical-magnetization vector representing concentration of molecules. We apply a convergence result from stochastic jump processes to deterministic trajectories and present a bifurcation analysis for the associated dynamical system. We show that non-symmetric module under study can exhibit very rich behaviours, including the empirical oscillations described by repressilator.


Introduction
Elowitz and Leibler [7] addressed the design and construction of a synthetic network providing the basic functionality of generating oscillations. Such functionality is essential to organize timemodulated biological functions like, for instance, the required periodic adjustment of an organism's physiology to the circadian rhythm [12,29]. Their idea was to construct a negative feedback loop composed by three transcriptional repressors that are not part of any natural biological clock. The resulting oscillating module was called repressilator (see Figure 1(a) for illustration).
We remark, as done by Elowitz and Leibler [7], that the repressilator displayed noisy behavior, this fact must be related with stochastic fluctuations of its components. Indeed, Elowitz et al. [8] verified that gene expression is inherently variable, or noisy, due to random fluctuations in individual cells (see also, Shahrezaei and Swain [31] and Swain et al. [33]).
Progress in the understanding of this module, as well as others naturally occurring networks, and their associated control mechanisms demand the development of mathematical models that manage good balance between simplicity and usefulness. In the literature there exist several different approaches to this aim. For instance, in the case of the repressilator we could find works such as Chen and Aihara [4], Chen et al. [5] and Wang et al. [37].
This kind of interaction modules are widely known as feedback loops, which carry out specific functions in a cell, decomposing realistic cellular control networks [1,32]. The study of feedbacks loops have received extensive attention in recent years [3,24,28]. Particularly, the development of general frameworks can be found in Along [2], Tyson and Novak [36], and especially in Wang et al. [38]. In general, we are interested in understanding for the dynamical behaviour of the concentration of molecules involved in the interactions.
In this work, we propose an alternative approach, which consists in the application of a class of interacting particle systems (see [22]). The ideas are taken from the type-dependent stochastic spin system, proposed by Fernández et al. [13]. Essentially, the main feature of this approach is taking account the inherent random variability in gene expression. The potential applications in general biological signaling networks are discussed in Fernández et al. [13].
As an illustration of our approach, we propose the study of a particular cyclic-interaction loop, that we call non-symmetric clock module, which is a simple extension of repressilator. Thus, to address this kind of qualitative approach, we propose the application of a mean-field type-dependent Ising model dynamics. Although the Ising model was initially studied to understand the physical phenomenon of ferromagnetism [16], nowadays this model represents a useful tool in different areas such as image processing, neural networks or earthquake dynamics [6,30,35], among others.
In our modeling scheme, the dynamics of the type-dependent Ising models is projected onto associated continuous-time jump processes, called density-profile processes, which are random walks mapping the macroscopic evolution of the particle systems.
Fernández et al. [13] showed that for arbitrary but fixed time intervals, in the limit of a very large number of particles (thermodynamics limit), the evolution of these jump processes converges to time dependent functions satisfying a correspondent deterministic dynamical system. We will include a simpler and straightforward proof of the convergence from the stochastic trajectories of the density-profile to deterministic paths ruled by non-linear differential equations. Our technique is based on the work of Ethier and Kurtz [11], who characterized a class of Markov jump processes called density dependent population processes (see also Kurtz [19]). Therefore, we use the convergence result to study the dynamical behaviours of our non-symmetric clock module. Particularly, we will analyse the influence of parameters in the behaviours of our model. The characterization of the evolution of the associated dynamical systems, that is a bifurcation analysis, completes our work.
We remark that the approach introduced by Fernández et al. [13] allows us to analyse general classes of feedback loops. For instance, we can enumerate the class of feedback loops studied in [21], which are simple enough to be analysed in a theoretical framework, but admiting very rich dynamical behaviours. Moreover, all the loops have been found in transcriptional regulatory networks (Leite and Wang [21]).
As we will see in section 3.2 the approach proposes a non-reversible stochastic microscopic dynamics based on interacting particle systems and statistical mechanics ideas. We used a mean-field version as a suitable approximation of other Ising-type interactions. In this sense, Mendonça and de Oliveira [23] discussed properties of the stationary measure of a TDSIM with near neighbours interaction in the study of repressilator. Further investigations could focus in other feedback loops, exploiting applications and to better understanding of the type-dependent stochastic dynamics.
The paper is organized in the following way. In section 2 we introduce the non-symmetric clock module, which will be used as an illustration. Section 3 includes the definition of type-dependent spin models, and particularly the dynamics of the mean-field type-dependent stochastic Ising model. We construct the associated density-profile process and prove convergence to deterministic trajectories in section 4. Section 5 includes a bifurcation analysis of the dynamical system associated to our nonsymmetric clock module. Finally, section 6 contains some conclusions and further investigations.
2 Feedback loops: the example of a non-symmetric clock module In this section we introduce a specific feedback module, which will be used to explain the step by step of our modeling scheme. However, we remark that the stochastic approach in Section 3 can be applied for dynamical studies of general feedback loops. We use the notion of network motifs to design the repressilator and our non-symmetric clock module. A feedback loop motif is a simple representation of a transcriptional regulation network [21,36]. That is a cycle in a directed graph whose vertices (also called nodes) can represent concentrations of proteins or genes (see Figure 1). In this sense, the edges can represent either positive or negative interactions. In other words, a positive (resp. negative) edge implies that the component in tail vertex activates (resp. represses) the transcription of the gene of the element in head vertex.
We shall consider a simple feedback loop that we call non-symmetric clock module, it is based on the repressilator proposed by Elowitz and Leibler [7]. The repressilator is a three transcriptional repressor systems that are not part of any natural biological clock, and were used to build an oscillating network in Escherichia coli. In other words, this is a negative feedback loop which provides the basic functionality of generating oscillations. Which is biologically required in several contexts like in cell-cycle and in the setting up of circadian cycle. Figure 1(a) shows the representation of the repressilator. We denote by A, B and C, the LacI protein, the TetR gene and the CI gene, respectively. Accordingly, the first repressor protein, LacI from E. coli, inhibits the transcription of the second repressor gene, tetR from the tetracyclineresistance transposon Tn10, whose protein product in turn inhibits the expression of a third gene, cI from l phage. Finally, CI inhibits lacI expression, completing the cycle. In other words, the rate of change of component A density at each time depends only on C density in an inhibitory manner: the density of A tends to decrease if the concentration of C is high. A similar dependence holds between C and B and between B and A.
In this work, we propose the study of a non-symmetric cyclic-interaction module, see Figure 1(b). We will refer to A, B and C as being three types of molecules. We focus on this simple feedback loop because it is large enough to admit complex dynamical behaviours. The model includes the parameter J that measures the strength of the interaction and a parameter δ ∈ [0, 1], which allows us to distribute the interactions between each pair of the three components (A, B, C). Note that the interactions between neighbour components could be non-symmetric, but the global interaction holds the invariance.
The particular cases when δ is equal 0 or 1, are similar to the repressilator above defined. Although, it is important to stress that for positive or negative values of J, we have inhibition or activation cycles, respectively. We are interested in a characterization of dynamical behaviour for this module, as a function of parameter J and δ.
Therefore, we describe a microscopic stochastic approach to derive an associated dynamical systems, which only seeks to incorporate the essential qualitative information about the biochemical interactions. This approach is based in the work of Fernández et al. [13], which borrows ideas from interacting particle systems [22], in such a way that spins represent the internal states of components of the feedback loop.

Modeling setup: microscopic type-dependent dynamics
In this section we describe the type-dependent stochastic spin systems proposed by Fernández et al. [13] to study signaling biological networks. We focus our explanation by studying the dynamical behaviours of the non-symmetric clock module exposed in previous section. In particular, we will define a type-dependent stochastic Ising model (TDSIM). That is, a microscopic model with Isingtype interactions and having dynamical evolution defined by the type-dependent dynamics proposed in Fernández et al. [13].

A family of interacting particle systems.
The type-dependent stochastic spin models are a family of stochastic spin-flip systems introduced in Fernández et al. [13]. This family extends the usual definition of particle systems [22], to allow asymmetric dependence of rates on the energy function, the Hamiltonian. As a consequence of this asymmetric dependence, a particularity of these models is its non-reversible stochastic dynamics.
Next, we explain these models. We need to introduce a set of spin types T , of cardinal k, representing genes or proteins (for instance: LacI, TetR and CI, in repressilator above). Also, a vertex set V of a simple finite graph of order N = |V|, denoting spatial positions available for each one of the types. We call the ordered pair (i, n) ∈ Λ = T × V a site. Moreover, the set of internal states for each i ∈ T , will be denoted by S i = {a 1 , . . . , a s i }. Hence, the spin system has site-space Λ and configuration space Ω Λ = i∈T S V i . For a configuration σ ∈ Ω Λ we denote by σ(i, n) the value of the spin at site (i, n) (see Figure 2).

Figure 2:
Lattice representation of the site-space Λ for a type-dependent spin system. Each vertical line represents a type i ∈ T . At horizontal lines we see the spatial positions V. The central spin shows the notation σ(i, n) ∈ S i , for a given configuration σ ∈ Ω Λ .
The continuous-time evolution of these models is governed by a non-reversible Glauber spin-flip stochastic dynamics. In other words, we have a stochastic evolution in continuous time, for which only one particle flips at each transition. Moreover, the dynamics is non-reversible with respect to the Gibbs measure (defined in (3.6)).
The corresponding rates are determined in terms of a function H Λ : Ω Λ → R, the Hamiltonian of the spin system. To define the Hamiltonian, we denote then, the Hamiltonian of these models is determined by a family of interaction matrices J n,l : C ×C → R, one for each pair of spatial positions n, l ∈ V. Due to the original applications in signaling biological networks, these matrices J n,l [·; ·] are not assumed to be symmetric. Particularly, we say that J n,l [(i, a); (j, b)] indicates the strength of the influence that a spin at a site (i, n) ∈ Λ in internal state a ∈ S i has upon a spin at (j, l) ∈ Λ that is in internal state b ∈ S j . As we illustrated in the repressilator module (see Figure 1(a)), type A components act upon type B components, while the reciprocal interaction does not occur. Then, it is reasonable that the quantity J n,l [(i, a); (j, b)] may be noticeably different that the one in the opposite direction J l,n [(j, b); (i, a)]. However, we are not assuming asymmetry with respect spatial positions, because the interaction is only associated to types, that is for n, l ∈ V, i, j ∈ T , a ∈ S i and b ∈ S j . Consequently, as usual in statistical mechanics, the Hamiltonian is then defined by for each configuration σ ∈ Ω Λ .

Type-dependent stochastic dynamics.
We define the stochastic dynamics proposed by Fernández et al. [13]. This dynamics allows only single spin-flips, more precisely, single-site internal-state transitions. In other words, we assume that at each time only one molecule could be produced or degraded. In statistical mechanics notation, we say that we have a Glauber type stochastic dynamics. Formally, given a configuration σ ∈ Ω Λ , a site (i, n) and an internal state a ∈ S i , we denote σ a That is, a configuration for which we fix the value a for the spin at site (i, n). Thus, the energy cost for the transition from σ a (i,n) to σ b (i,n) is given by Usually, this total change of the energy associated to the flip at site (i, n) from state a to b, is used to define the rates of transition in particle systems. Hence, generally we have a continuous-time stochastic evolution being reversible with respect to the Gibbs measure, for each σ ∈ Ω Λ . However, in type-dependent spin models the definition of the rates of transition is quite different. The asymmetry of the interaction defined above leads naturally to the decomposition of (3.5) in the following manner, where, that collects the change in the influence of the configuration σ upon the site (i, n), when internal state there changes from a to b. On the other hand, collects the change of the influence that the site (i, n) has on all other sites when its internal state flips from a to b. Furthermore, as a particularity of the type-dependent spin models, each transition rate depends only on the energy changes brought upon the site (3.8). Thus, we denote λ a→b (i,n) (σ) the rate of a transition flipping σ a (i,n) to σ b (i,n) , which depends only on (3.8), that is, Finally, we state a formal definition of the spin models proposed by Fernández et al. [13].
on Ω Λ , defined by a spin model with a type-dependent interaction, given by the Hamiltonian in (3.3) and a dynamics with rates of the form (3.10).
The dynamics of these spin systems yields non-reversibility with respect to the Gibbs measure (3.6). Then, the type-dependent dynamics has mathematical interest for itself. It is natural to study the stationary measure for particular cases of the interactions matrices (J i,j ) i,j∈C . That is, for instance local interactions like nearest neighbors (n.n), as usual in many particles systems. In Mendonça and de Oliveira [23], was proposed an Ising-type n.n interaction to analyse, by Monte Carlo simulations, some properties of the stationary measure for the repressilator.
In this work, as done in Fernández et al. [13], we consider a mean-field interaction for an Ising model. The mean-field interaction is a natural approximation of local interactions and given our interest in concentrations of molecule, as global observable, this is a suitable choice. Of course, the convergence results to be explained in Section 4 hold whenever we assume mean-field dynamics.
where {α i,j } i,j∈C , is a real matrix.
We conclude this section explaining the particular modeling setup that we will use to study the non-symmetric clock module.

Mean-field TDSIM
We follow the ideas of Fernández et al. [13] to model our non-symmetric clock module through a type-dependent stochastic Ising model (TDSIM), which has all internal state spaces S i = {−1, +1}, for all i ∈ T .
We think the spin types as points {A, B, C} over the circle (see Figure 1(b)) and, for each i ∈ T = {A, B, C}, we denote h(i) the neighbour of i in the clockwise direction. And a(i) to be the neighbour in the anti-clockwise direction. Borrowing statistical mechanical nomenclature, we say that for J < 0, the activation interactions in the Figure 1(b) are ferromagnetic. Moreover, for J > 0, the inhibition cycle is said to be an anti-ferromagnetic model.
Next, observe that the set (3.1) will be particularly defined by In addition, the set of vertices V will represent a reservoir of capacity N (actually, we have three reservoirs, one for each type of molecule in T ). In this way, for a configuration σ ∈ {−1, +1} Λ , the spin value σ(i, n) = +1 means that there is a molecule of type i ∈ T at spatial position n ∈ V. Otherwise, if σ(i, n) = −1, it means the absence of molecules of type i ∈ T in the corresponding reservoir at spatial position n ∈ V.
Acordingly, the mean-field TDSIM must be defined by rates functions as in (3.10), where Φ(∆) = e −∆ (for details, see [13]). In addition, the real matrix {α i,j } i,j∈C * as in (3.11) is given by for (i, a), (j, b) ∈ C * , and where J ∈ R, δ ∈ [0, 1], and κ i ∈ R is an external field (or chemical potential). It must be clear that a site (i, n) ∈ Λ in internal state −1, will not influence others sites, because at (i, n) there is no molecule.
Notice that, as we said above, for δ = 0 (as well as δ = 1) and J > 0, the totally asymmetric case, we obtain an Ising model for the repressilator represented in Figure 1(a).
Finally, the transition rates for TDSIM are defined for each σ ∈ Ω Λ by where a + i = |{l ∈ V : σ(a(i), l) = +1}|, and h + i = |{l ∈ V : σ(h(i), l) = +1}|, for i ∈ {A, B, C}. For each type i, the parameter κ i must be interpreted as the own rate of production of molecules i ∈ T .
In general κ i is a linear function of parameter J. Because we are interested in the relation between two-body interaction (J) with individual interaction (κ).
We summarize our modeling setup in the following statement: The microscopic evolution of the non-symmetric clock module is described by a mean-field type-dependent stochastic Ising model with continuous-time Glauber dynamics (σ t ) t≥0 , whose rates of transition are given by (3.14).
In the next section we will focus in the macroscopic evolution of the non-symmetric clock module. Thus, we shall define the density-profile processes, which are a family of random walks with continuous time evolution. Moreover, we will prove an almost sure convergence from these stochastic trajectories to deterministic paths governed by non-linear differential equations.
4 Macroscopic evolution: the density-profile process and its convergence to deterministic dynamics In this section we focus our attention in the concentrations of biochemical components into our feedback loop defined in Section 2. Hence, we will mostly be interested on the vector of empirical magnetization of the mean-field TDSIM. Finally, we will study the thermodynamic limit, that is, the behaviours when |V| = N → ∞.
We project the mean-field TDSIM onto a jump Markov process in R 3 , called density-profile process, that represents the densities of the three biochemical components in our feedback loop. Thus, in statistical mechanics notation, the stochastic spin system and the density-profile process provide, respectively, the microscopic and macroscopic views of the same model. where We describe the evolution of that processes in the following way: for each position X N (t) = x ∈ D N , the jumps of the form where l ∈ J , are defined by a collection β l (x) of functions, β l : D N → [0, ∞], l ∈ J . Of course, we require that x ∈ D N and β l (x) > 0, imply x + l/N ∈ D N . It must be clear that each β l (x) will be defined as a function of the transition rates in (3.14). Therefore, the jump in (4.3) occurs with rate where x i = X N i (t) as defined in (4.2), and e i is the unitary vector in the direction of i ∈ {A, B, C}. Note that each variable x i represents the density of elements of type i, or we say as well, the concentration of the component i ∈ {A, B, C} in our non-symmetric clock module. Thus, N x i will represent the number of molecules of the type i, and in some abstract sense, N (1 − x i ) denotes the available number of possible molecules into the reservoir of size N . Although, in the Ising system, this last quantity is just the number of sites of type i ∈ T , with spin-value equal to −1.
The next step is the characterization of these jump processes in the thermodynamics limit, that is N → ∞. Particularly, Fernández et al. [13], showed that in the limit of a very large number of spins, these density-profile processes converge to time dependent functions satisfying a correspondent deterministic dynamical system. They used a pathwise approach, which strongly exploits large deviation theory [26], and coupling of random variables [34]. Their method provides explicit bounds for the distance between the stochastic and deterministic trajectories.
However, the qualitative study of the behaviors of our feedback loop, does not require this kind of accurate results obtained in [13]. Therefore, we follow the works of Kurtz [19] and Ethier and Kurtz [11], to obtain a simpler and straightforward proof of the convergence from the paths of the density-profile processes to deterministic trajectories, these latter ruled by non-linear differential equations. Now, we aim to state our first theorem and prove it. The main tool that we will use is the characterization of a family of jump processes given by Kurtz [19]. We stress that our result can be easily extended to more general models, that include applications to a very extensive class of feedback loops.
Initially, we define a jump Markov process X N on Z 3 , for which we will use some properties given by Kurtz [19]. Particularly, the jumps of this process are of the form l ∈ J (see Definition 4.1), and their intensities will be given by β * l (k) = N β l (k/N ), k ∈ Z 3 , where β l (x) as defined in (4.5). Thus, let X N (0) being non-random, the evolution of X N at time t ≥ 0, can be written as where Z l (t) := |{s ≤ t; X N (s) − X N (s−) = l}|, that is, the number of jumps of size l until time t. Then, by some results from Chapter 7 of [19], we have: for each t ≥ 0, where the Y l (u), l ∈ J , represent independent Poisson processes with corresponding intensities u. Now, setting (4.8) and X N = N −1 X N . Thus, we have is the Poisson process centred at its expectation. Thus, our first result is the almost sure convergence of stochastic trajectories (4.1) to associated deterministic paths. The proof of this result is essentially repeat the arguments of Theorem 2.1 from Chapter 11 of [11], but we include the proof for completeness.
Proof. First of all, denoteβ l := sup x∈D N β l (x). Thus, it is easy to see that l |l|β l < ∞, (4.12) and by the Lipschitz continuity of the rates of the jumps (4.5), there exists a fixed M > 0 such that Observe that by (4.9) and (4.10), we have for each t ≥ 0, (4.14) Now, let denote Therefore, using (4.13), (4.15) and the Gronwall inequality, we can rewrite (4.14) as follows, for all t ≥ 0. Thus, since ε N (t) is a non-decreasing function, To obtain (4.11) we need the following lemma. Proof. First, note that we have the following equivalence in distribution where Y i l (u) are independent Poisson processes with rate u. Then, the strong law of large number implies that Furthermore, by definition ofβ l above, it follows that If u ≤ t, a basic property of Poisson processes states that, for each l ∈ J . Then, 23) for all N ≥ 1 and each l ∈ J . Now, note that, where Y i l (u) are again independent Poisson processes with instensity u. Finally, by the law of large numbers applied to the independent random variables in bracket at r.h.s of (4.24), and by condition (4.12), That is, by (4.23) we can interchange the limit and summation for the expression at r.h.s in (4.21). Therefore, using (4.20), the Lemma follows. Next section is devoted to our results about modeling a specific biological feedback loop, that is, our non-symmetric clock module. Clearly, the application of the Theorem 4.2 allows us to analyse the qualitative behaviour of the concentrations of the molecules involved in the interactions. Therefore, we are able to include the role of stochasticity in the gene expression, but we also simplify the qualitative study through a bifurcation analysis for the associated dynamical system.

The associated dynamical system: formulation and bifurcation analysis
In this section we include a qualitative study of the dynamical behaviours of the non-symmetric clock module. We use the ideas of previous sections to analyse a dynamical system related to stochastic evolution of the concentration of molecules involved in our feedback loop. First, we briefly summarize the ideas in previous sections. Initially, we stated a mean-field TDSIM (σ t ) t≥0 on Ω Λ and defined its dynamics by the rates of transition (3.14). Thus, the concentrations of each component in our non-symmetric loop (types A, B and C) were characterized by a densityprofile process, defined in (4.1). In the Ising case, the only independent variables are the densities of activated types that we denoted X N i (t), i ∈ {A, B, C}. The rate of a jump in density-profile process was defined in (4.4) and (4.5), based on the rate of corresponding transition in TDSIM. After that, in Section 4, we also studied the thermodynamic limit of the particle system (N → ∞), proving that the density-profile process converges to deterministic trajectories governed by non-linear differential equations (see Theorem 4.2).
Finally, in this section we show that depending on the parameter values, the magnetization random path can either converges to a unique stable fixed point (if −1 < J < 2), converges to one of a pair of stable fixed points (for J < −1), or asymptotically evolves close to a deterministic orbit in R 3 (when J > 2). Therefore, we follow Theorem 4.2 to study the dynamics of our non-symmetric clock module by analysing the associated dynamical system. However, we remark that the behaviours of the dynamical system are deterministic approximations of the stochastic evolutions of the jump processes in our approach.
The system of differential equations associated to the cycle-interaction module in Figure 1(b), and with rates given by (3.14) will be expressed bẏ for i ∈ {A, B, C}, J ∈ R, δ ∈ [0, 1], and κ i ∈ R. As usual, we will consider κ i = J/2, i ∈ {A, B, C}, that it, proportional to the molecules interaction parameter. The particular value is to reduce the number of parameters and to obtain suitable fixed points in dynamical analysis (steady state of concentrations (1/2, 1/2, 1/2)). Our next result is the statement of a bifurcation analysis to guarantee for the non-symmetric clock module that: for 0 < J < J c = 2, the concentration of all three components (A, B, C) remains stable, but oscillates with large amplitude as soon as J increases past the threshold J c = 2. For J < 0, we can see the appearance of two stable points when J < −1. These situations are showed in Figure 3.  Proof. First, it is easy to see that the dynamical system has a fixed point at x 0 = (1/2, 1/2, 1/2) T , becauseẋ = F (x 0 , J, δ) = (0, 0, 0) T , for all J ∈ R and δ ∈ D.
From now on, we consider the dynamical systemẏ = F * (y, J, δ), where y = (y 1 , y 2 , y 3 ) T and x i = 1 2 + y i . Then, the fixed point for the new system is 0 = (0, 0, 0) T . Following the theory in Perko [27] and Kuznetsov [20], we could guarantee that near 0 = (0, 0, 0) T , the dynamical system is locally topologically equivalent to its linearisationẏ = Ay, where A denotes the Jacobian matrix dF * dy evaluated at 0. Then, y 0 is stable if all eigenvalues λ of A satisfy Re λ < 0, and stability is lost when one of the real parts becomes positive.
Notice that A = ∂F * (0) is given by the expression The eigenvalues of this matrix are , and Therefore, the fixed point y 0 loses the stability at J c = −1, when λ 1 crosses the imaginary axis through the origin. Thus, we have two stable fixed points. This bifurcation is called fold (or pitchfork) bifurcation. The values of these fixed points could be calculated by solving: sinh(2Jy) + 2y cosh(2Jy) = 0, thus we obtain the relation J = 1 4y · log 1 − 2y 1 + 2y (5.4) for y = 0. We notice the symmetry around J-axis, this fact could be seen at left hand of Figure 3.
On the other hand, the stability is lost when λ 2 and λ 3 , symmetric around the real axis, cross the imaginary axis. This fact occurs when J = 2. The phenomenon includes the appearance of a stable orbit, this kind of bifurcation is known as Hopf bifurcation.
The oscillations experimentally verified in the repressilator (Elowitz and Leibler [7]) are explained in this dynamical behaviour through the (supercritical) Hopf bifurcation with respect to the interaction strength parameter J.

Deterministic macroscopic evolution in details
Let us finally show a detailed analysis of the dynamical behaviours of the system using linearisation tools. In fact, we do not study the dynamical system (5.1), however we analyse the linear system that approximates it.
We will exhibit a qualitative analysis of the concentration of the molecules involved in the interactions. The idea is a better understanding of the influence of parameter δ in the dynamical behaviours.
Equivalently, remember that y = (y 1 , y 2 , y 3 ) T and x i = 1 2 + y i , for i = {1, 2, 3}. We can say that the system Z 1 Z 2 Z 3 is obtained by the following two operations in the system Y 1 Y 2 Y 3 : Finally, to obtain the new coordinates of the system, we have (z 1 , z 2 , z 3 ) = (y 1 , y 2 , y 3 ) · R, with Again, by the theory in [20,27], the behaviour of the dynamical system in Z 1 Z 2 Z 3 , could be studied by linearisation, that is, we analyse where A = R −1 · ∂F * (0) · R, because R is an orthonormal matrix. Then, and expression (5.6) represents the following system: Notice that the behaviour over Z 3 -axis is independent of the dynamics on Z 1 Z 2 -plane. Then, we could establish the following two elementary conclusions: C1) For any J > −1, the trajectory on Z 3 -axis goes to zero, becauseż 3 is negative for z 3 > 0, and positive for z 3 < 0; C2) For any J < 2, the dynamics on Z 1 Z 2 -plane goes to zero in both directions. In other words, it goes to the point (0,0). Of course, the velocity and direction depend on δ and J.

Conclusions and comments
On one hand, based on a particular feedback loop we have described a theoretical framework to model and analyse the non-linear dynamics of gene regulatory networks. In the literature there exist several different approaches to this aim. However, in previous sections we showed that our approach, based on the well studied Ising model, can include the inherent variability in gene expression, also this microscopic approach could capture detailed characteristics of the system interpreted as a promotioninhibition circuitry. Thus, it may be particularly useful to illuminate how simple feedback loops manage to perform their basic functionalities and how these functionalities are further integrated into the whole cellular regulatory network.
On the other hand, the mathematical study of the type-dependent dynamics suggests a longstanding theme of statistical physics of nonequilibrium systems. That is, the question of the nature of their stationary measure, and in particular of the existence of phase transitions. Thus, it will be interesting for us to propose local interactions in type-dependent Ising models, instead of mean-field interaction, at this point we could consider several lattices (like Z d or triangular lattice) and treat these questions.
Particularly, the study of the TDSIM with local interactions is closely related with works like Godrèche [14] and Godrèche and Bray [15] (see also de Oliveira [25]). They considered a kind of asymmetric Ising dynamics, called directed Ising model. In [14,15] was studied several dynamics, like Glauber or Metropolis. They obtained a large space of parameters defining the rates functions allowing irreversible Gibbsian Ising models, whenever the dynamics is not totally asymmetric.
Furthermore, bifurcations showed in the associated dynamical systems (Section 5) suggest the existence of phase transitions in the type-dependent Ising model. Of course, it would also be natural to consider metastability issues [26] to characterize the dynamical behaviors of this particle system. In this sense, a related work is due to Kotecký and Olivieri [17], they studied a discrete-time Metropolis dynamics for a two dimensional ferromagnetic asymmetric Ising model, in which the vertical and horizontal interaction parameters have different values.
Finally, we say that another very interesting mathematical question is the characterization of loss and recovery of Gibbsianness in stochastic evolutions, considered in recent works, in the area of mathematical physics. For instance, van Enter et al. [9] studied some Ising models for which this phenomenon occurs. In addition, Kulske and Le Ny [18] initiated a fruitful research direction showing that Gibbs-non-Gibbs transitions can also be defined naturally for mean-field models.