On Hybrid Gene Regulatory Networks
Adrian Wurm Honglu Sun








A. Wurm

BTU Cottbus-Senftenberg, Lehrstuhl Theoretische Informatik,
 Platz der Deutschen Einheit 1, 03046 Cottbus, Germany,

<https://www.b-tu.de/>
wurm@b-tu.de Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France



On Hybrid Gene Regulatory Networks
    Adrian Wurm^1 Honglu Sun^2 
    April 29, 2024
==================================



    In this work, we study a class of hybrid dynamical systems called hybrid gene regulatory networks (HGRNs) which was proposed to model gene regulatory networks. 
    In HGRNs, there exist well-behaved trajectories that reach a fixed point or converge to a limit cycle, as well as chaotic trajectories that behave non-periodic or indeterministic. 
    In our work, we investigate these irregular behaviors of HGRNs and present theoretical results about the decidability of the reachability problem, the probability of indeterministic behavior of HGRNs, and chaos especially in 2-dimensional HGRNs.


 



§ INTRODUCTION

Hybrid systems are dynamical systems that admit both continuous and discrete behavior. 


They were initially proposed to model the interactions between the discrete digital world and the continuous physical world <cit.>.
In recent years, hybrid systems were also used in system biology to study biological networks<cit.>.

In this work, we study a class of hybrid systems called hybrid gene regulatory networks (HGRN) <cit.>, which is an extension of Thomas' discrete modeling framework <cit.>. 
These hybrid systems were proposed to model gene regulatory networks, a framework to describe the regulation influences between genes.

Two widely used formalisms to study gene regulatory networks are discrete models (like Boolean networks <cit.>) and continuous models (differential equations <cit.> and stochastic models <cit.>).
The major advantage of HGRNs compared to discrete models is that they can provide and process temporal information of the systems.
Compared to purely continuous models, the dynamical properties of HGRNs are easier to analyze, which makes them a more reasonable choice to apply on large and complex instances.

HGRNs are similar to piecewise-constant derivative systems (PCD systems) <cit.> which are a special case of hybrid automata <cit.>. The major difference between HGRNs and PCD systems, as for example in <cit.>, is the existence of sliding modes, which means that when a trajectory reaches a black wall (a boundary of the discrete state which can be reached but cannot be crossed by trajectories), the trajectory does not end there but is instead forced to move along this black wall.
There exist other methods to define the behavior of trajectories on a black wall <cit.> which are different from the sliding modes in HGRNs.

One important research topic of dynamical models of gene regulatory networks is how to analyze the reachability of states, for instance, whether a trajectory starting in a certain state can reach another given state.
These questions arise in cell reprogramming <cit.>. In <cit.>, the reachability problem of HGRNs is investigated and methods are proposed to analyze the reachability of regular trajectories, meaning those that reach a fixed point or converge to a limit cycle. 

In HGRNs however, there also exists so-called chaotic trajectories that cross an infinite sequence of discrete regions so that this sequence does not end up in a cycle.
For now, our understanding of this irregular behaviour is quite limited.
Further understanding of these trajectories can potentially lead to the development of new methods to analyze the reachability problem for them. 
In this work, we want to fill the gaps in this research field by investigating irregular behavior occuring in HGRNs. 
Our major contributions are:

    
  * We prove that the reachability problem for HGRNs is undecidable and that this enforces the existence of chaotic trajectories.
    
  * We prove that under reasonable conditions, the probability for a random trajectory to admit indeterminism tends to 1 as the dimension of the system tends to infinity.
    
  * We argue that  in two-dimensional HGRNs, chaos is not possible without indeterminism, and conclude that the probability for a random two-dimensional HGRN to admit chaos with indeterminism is 0. 


The paper is organized as follows:
In Section <ref> we collect basic notions and recall the definition of HGRNs as well as their properties. In Section <ref> we prove that HGRNs are Turing-powerful and conclude that chaotic trajectories must exist. In Section <ref> we give some structural results on indeterminism and provide conditions as well as probability estimations for the bifurcation of trajectories.

The paper ends with some open questions. 







§ PRELIMINARIES AND NETWORK REACHABILITY PROBLEMS



We start by defining the problems we are interested in; here,
we follow the definitions and notions of <cit.> for everything related to 
hybrid gene regulatory networks. 

Intuitively, a HGRN is a  an N-dimensional rectangular cuboid/hyperrectangle with integer edge lengths decomposed into hypercubes with unit edge length. Each such hypercube d is called a discrete state and is associated with a celerity c_d enforcing the direction in which a trajectory is allowed to move inside d, so a trajectory may only change its direction in the moment of transition between two discrete states.



Let N∈ℕ and let n_i∈ℕ for each i∈{1,...,N}. A Hybrid gene regulatory network (HGRN) on N genes with n_i levels is a pair ℋ=(E, c), where

    E := { d ∈ℕ^N |∀ i ∈{1,2,...,N } : 0≤ d_i ≤ n_i }




is called the set of discrete states, and c:E→ℝ^N is called the celerity function.



A hybrid state of a HGRN is a pair h=(π,d) containing a fractional part π∈[0,1]^N, and a discrete state d∈ E.
E_h is the set of all hybrid states and E_sh is the set of all finite or infinite sequences of hybrid states: 

    E_sh = { (h_0,h_1,...,h_m) ∈ (E_h)^m+1| m ∈ℕ∪{∞}}.

Note that E_h=:E_h^1⊆ E_sh.



A boundary of a discrete state d is a set of hybrid states defined by e(i,π_0,d) := { (π,d) ∈ E_h |π^i = π_0, }, where i ∈{ 1, 2,..., N} and π_0 ∈{0,1}.



 Let t_0 ∈ℝ^+ ∪{∞} and [0,t_0] a time interval. For a function τ:[0,t_0]→ E_sh, we denote by abuse of notation by π:[0,t_0]→ E_sh the projection of τ on the second component. 
 
 
 For τ(t)=h, we denote by  dτ(t^+)/dt:=dπ(t^+)/dt the temporal derivative of the fractional part of τ, dτ(t^-)/dt and dτ(t)/dt are defined analogous.
 
 We call τ a hybrid trajectory, if for every time t∈[0,t_0] with τ(t)=h and h=(h_0,...,h_m)∈ E_h^m+1 a sequence of m+1 hybrid states h_j=(π_j,d_j), the following conditions are satisfied:




  * For t=0 we have m=0, and τ(0)∈ E_h is not on any boundary.

  * If h_0 does not belong to any boundary, then m=0 and dτ(t)/dt = c(d_0) for t>0 and dτ(t^+)/dt = c(d_0) for t=0 

  * If h_0 only belongs to one boundary e, let e be the upper boundary in the i^th dimension (the following is analogous when e is the lower boundary). In case d_0^i is not the maximal discrete level of gene i, let the discrete state on the other side of e be noted as d_r, where d_0^k = d_r^k for all k ≠ i, and d_0^i+1 = d_r^i.
There are four possible cases:
	

	
  * If c(d_0)^i = 0, then m=0 and the trajectory from the current state will slide along the boundary e with speed dτ(t^+)/dt = dτ(t^-)/dt = c_0, which is then called a neutral boundary of d_0. 
	
  * If c(d_0)^i > 0, and either d_0^i is the maximal discrete level of the i^th gene, or d_0^i is not the maximal discrete level of the i^th gene but the i^th component of c(d_r) is negative, then m=0 and the trajectory from the current state will slide along the boundary e, which is called an attractive boundary of d_0. We then have that dτ(t)/dt^k = c(d_0)^k for all k ≠ i, dτ(t^+)/dt^i = 0 and dτ(t^-)/dt^i ∈{0,c(d_0)^i}.
	
  * If c(d_0)^i > 0, d_0^i is not the maximal discrete level of the i^th gene, and the i^th component of c_r is positive, then the trajectory from the current state instantly crosses  the boundary e and enter the discrete state d_r. We then have m=1 and τ(t)=((π_0,d_0),(π_1,d_1)), where d_1=d_r, π_0^j=π_1^j for all if j ≠ i, and π_1^i = 0. Such a boundary is called an output boundary and the counterpart of d_0 is called  an input boundary of d_r. We then demand that dτ(t^+)/dt = c(d_r) and  dτ(t^-)/dt = c(d_0).
	

  * If h belongs to several boundaries, then the previous cases can be mixed:
	
	
  * If in these boundaries there is no output boundary, then the trajectory from the current state will exit all input boundaries and slide along all attractive or neutral boundaries.
	
  * If in these boundaries there is only one output boundary, then the trajectory from the current state will cross this output boundary.
	
  * If in these boundaries there are several output boundaries, then the trajectory from the current state can cross any of them, but can only cross one boundary at one time, which causes indeterministic behavior.
	


A hybrid trajectory τ is called 
    
  * finite, if t_0≠∞ and the sequence of states that are reached is finite,
	
  * periodic, if it is defined on t_0=∞ and ∃ T, t_1>0, ∀ t ∈ [t_1, ∞[, τ(t) = τ(t+T), 
    
  * chaotic, if the sequence of discrete states that it passes is infinite and non-periodic, meaning it does not end up in a cycle and
    
  * regular, if it is not chaotic. 














 












The reachability problem for HGRNs is the following: Given a HGRN and two sets S,T of hybrid states in form of linear inequalities, the question is whether there exists a trajectory τ starting in a hybrid state s∈ S and reaching at some time t a state f∈ F, meaning 

    ∃τ∃ s∈ S,f∈ F, t>0 : τ(0)=s τ(t)=f.

If such a trajectory exists, we say that F is reachable from T.

Note the special cases that S,T might contain only one hybrid state or that they may be a whole discrete state. It was shown in <cit.> that the reachability problem is decidable for regular trajectories.


Consider a subset A of hybrid states of a HGRN. The area of attraction φ(A) is the set of all hybrid states with the property that every trajectory starting in them will eventually intersect A. The extended area of attraction Φ(A) is the set of all hybrid states with the property that there exists a trajectory starting in them that will at some point intersect A. 



 Fig. <ref> is a sketch of a HGRN on three genes with levels n_1=2,n_2=4,n_3=2, the arrows indicate the celerities where an encircled circle indicates an upward movement and an encircled X a downward movement, missing arrows indicate that the celerity is the zero-vector. The left side of the figure illustrates level 1 of gene 3, the right side level 2.
 

 The network computes the binary expansion of a number in [0,1] in the following way: If a trajectory passes A=e(2,1,(1,2,1)) (blue area), interpret its distance to the rightmost part A∩ B , B=e(1,1,(1,2,1)) as its current value x. If x<1/2, the trajectory will enter the small circle of discrete states (1,2,1),(2,2,1),(2,1,1) and (1,1,1) and finally end up in A again with value 2x. If x>1/2, it will pass the big circle of  discrete states (1,2,1) (where 1/2 is subtracted from the value),(1,3,1),(1,3,2),(1,4,2) (where its value is multiplied with two) ,(2,4,2),(2,3,2),(2,2,2),(2,1,1) and (1,1,1) and again reach A, this time with value 2(x-1/2). If one follows the trajectory and interprets the first case as 0 and the second case as 1, one obtains the binary expansion of the initial x.

 If for example we start in a=((3/16,3/16),(1,1,1)), we obtain the trajectory drawn in red. It passes A at value 3/8, so it runs through the small circle and enters A again at value 3/4, therefore entering the big circle ending up in A again at value 1/2. So far, we obtained the binary expansion 0.01, but now the trajectory becomes indeterministic because the celerities allow a bifurcation at ((1,1),(1,2,1)). This corresponds to the fact that 3/8 has the two binary expansions 0.011 and 0.0101̅.

 Observe that we obtain a periodic trajectory for every rational value x and a chaotic trajectory for x irrational.




§ TURING COMPLETENESS AND CHAOS




Turing-completeness was already shown for PCD systems in <cit.>, a model related to HGRNs. We would now like to prove it for HGRNs, for they are more restrictive in the sense that they only allow a change of derivatives at integer coordinate changes and admit indeterminism with positive probability. 


HGRNs are Turing powerful.




Let M=(Q,Σ,Γ,δ,q_1,□,F) be a deterministic Turing-Machine starting on an empty tape, Q={q_1,...,q_n} , F={q_n}. Assume that all symbols are encoded in binary, so Σ=Γ={0,1} and that the encoding of □ consists only of zeroes. Also assume that for each two states q_i and q_j there are transitions in only one direction, either q_i→ q_j or q_j→ q_i but not both, a property that can be guaranteed by introducing auxiliary states. We construct an equivalent HGRN H with the following genes:



  ∙ One with two discrete levels storing the current value called G_c 

  ∙ G_ℓ,1 with three discrete levels and G_ℓ,2 with two discrete levels for the tape values left of the head, where the current value is stored as the binary expansion with the most significant bits corresponding to the values closest to the head in the fractional part of the upmost discrete state of G_ℓ,1.

  ∙ G_r,1 and G_r,2 analogous for the tape values right of the head.

  ∙ One called G_i for every q_i∈ Q with two discrete levels each.

  ∙ Two called G_α and G_β with two discrete levels each, keeping track of the current calculation step as in Fig. <ref>.

The whole network is therefore of shape [0,2]× [[0,3]×[0,2]]^2×[0,2]^n×[0,2]^2 and our trajectory starts in ((0,0),[(2,0),(0,0)]^2,(1,1),[(0,0)]^n-1,(0,1),(0,0)).

Let (q_i,a)→(q_j,b,d) be a transition in δ, without loss of generality d=r stating that the head moves to the right (d=ℓ is analogous and d=n is easily derived). We choose the celerities as follows:





In the first calculation step α=β=0 we put q_j into the upmost position. Note that we do not loose the information on which state we started in because we assumed that there is no transition from q_j to q_i.




The celerities describing how to update the left side tape start with making the bits less significant by one position because the head recedes by one position, which corresponds to dividing the numerical value of G_ℓ,1 by two, see Fig <ref>, left. This is done in the first calculation step α=β=0 along with adding half of b to it, see Fig. <ref>, middle, because b is on the position left of the head after that transition.



In the second calculation step α=0,β=1, this value is then copied back from G_ℓ,2 to the original position, see Fig. <ref>, right, G_ℓ,1, both G_ℓ,1 and G_ℓ,2 remain inactive for the last two steps.



The celerities describing how to update the right side tape work similar. We start by multiplying the value by two, see Fig. <ref> left, because each bit now becomes more significant by one position. This is done in the first calculation step α=β=0, where we also catch the trajectory on the border between G_r,1=1 and G_r,1=2 if it is in G_r,2=1 and let it pass until the lower boundary of G_r,1=1 if G_r,1=2, in order to keep track of what our new current value G_c will be in the next Turing step. 




Note that if the initial value was 1/2, we have an indeterminism when hitting the lower boundary of G_r,1=3 and the upper boundary of G_r,2=1. In the second calculation step α=0,β=1, we therefore adjust the celerities to make sure that we end in G_r,2=2 whenever the initial value was greater than or equal 1/2, see Fig. <ref> middle (the case that all positions on the right side tape are 1's is not possible).

In the third and fourth calculation step α=1,β=1 and α=1,β=0, we pop the most significant bit of the right side tape to be the new G_c, which means that if the value of G_r,1 is greater than 1 we subtract 1 from G_r,1 and add move G_c to the upmost position, otherwise to the lowest, see Fig. <ref> right. 





Note that we loose the information on whether we moved left or right, so we actually have the celerities of the fourth step on the right side of the tape for the left side of the tape as well. However, the left side is not effected by that because it stays in  G_ℓ,2=3 from the end of step two on and therefore never enters a hypercube with non-zero celerities induced be the update of G_c. 

In the fourth step we also bring G_i to the lowest possible point. By assumption of anti-symmetry, i can still be distinguished from j at that point without the knowledge of d.▪




The reachability problem for HGRNs is undecidable.



The halting problem can be reformulated as the reachability problem of state q_n.▪


Note that the trajectories that are obtained by the Turing machines that behave chaotic automatically lead to chaotic trajectories, for the sequence of reached states cannot have a repeating pattern.  



§ INDETERMINISM

Most attempts to solve verification tasks such as the reachability problem become a lot more difficult as soon as trajectories are not regular any more. It is therefore desirable to estimate how likely a trajectory is to be chaotic or indeterministic.

In opposition to what is claimed in <cit.>, indeterminism will happen with non-zero probability in networks with at least 3 genes having at least two discrete levels each when all celerities and the starting point are chosen uniformly random in (0,1). Fig. <ref> is a sketch of a HGRN with 3 Genes G_1,G_2 and G_3 with two levels 0 and 1 each that allows the indeterministic red trajectory starting in point a in (0,0,0).




The trajectory reaches the attractive boundary to discrete state 1 of G_1, slides along until it reaches the attractive boundary to discrete state 1 of G_2 and slides "upwards" along this intersection of discrete domains until it reaches the hybrid state ((1,1,1),(0,0,0)). The border to the discrete state 1 of G_3 is an output boundary, so it is crossed and the hybrid state ((1,1,0),(0,0,1)) is immediately reached. Now the boundaries to the discrete levels 1 of G_1 and G_2 both are output boundaries, so the trajectory may pass any of them, leading to indeterminism. Note that we did not rely on probability-zero-events; The probability of the celerity in the first hypercube to point into the middle in all three dimensions is (1/2)^3, the first two reached boundaries are then attractive with probability 1/2 each as well as the third one being and output boundary and the two new boundaries after the crossing are output boundaries with probability 1/4 each. The accumulated probability is then 1/8·(1/2)^2·1/2·1/4^2=(1/2)^10. This probability significantly rises once we increase the size of the network, meaning the number of genes as well as their dimensions.


We can be more general:



Let all celerities of a HGRN with n genes and at least two discrete levels per gene as well as the starting point of a trajectory be randomly chosen by the uniform distribution in [0,1]^n.
Then the probability for this trajectory to lead to indeterminism tends to 1 as n tends to ∞. 




First, note that for each gene we have that an inner boundary is reached at some point with probability at least 1/2. These events are independent, so it follows by the central limit theorem that the probability of at least one fourth of all genes reaching an inner boundary tends to 1. Each of these n/4 boundaries is attractive with probability 1/2, if we group them into n/12 triples of boundaries, each triple starts with two attractive boundaries followed by an output boundary with probability 1/8 as in Example <ref>. The probability that after the crossing of such a third boundary, the two other boundaries become output boundaries is (1/2)^2 each, so the accumulated probability for each triple to lead to indeterminism is at least 1/8·(1/2)^2·(1/2)^2=(1/2)^7. The probability that none of the n/12 triples leads to indeterminism is therefore at most (2^7-1/2^7)^n/12 0.▪


Note that Example <ref> relied on the third dimension, for hitting two edges at the same time without sliding before is highly unlikely:


Random HGRNs in two dimensions admit chaos only with probability 0.




Since a trajectory cannot cross itself, it is obvious that chaos can only be obtained by bifurcation.

A trajectory starts in the inner of a hypercube with probability 1. In order to bifurcate, it needs to reach two output boundaries at the same time, meaning it has to hit one of finitely many points. Reaching such an intersection point q by sliding towards it along an attractive boundary will not lead to indeterminism at q because the boundary which one is sliding along is not an output boundary. A point p in a two-dimensional HGRN that has an area of attraction with non-zero measure must either have integer coordinates connected to an attractive boundary or be reachable from such a point, which is a zero-set.

The area of attraction of a point that admits two output boundaries is therefore one-dimensional with probability 1, meaning it is a zero-set.  ▪





Note that an indeterministic network such as in Example <ref> inherently enables chaotic trajectories if both choices lead back into the area of attraction of the point a at which the path was bifurcating, therefore forming two cycles T_1,T_2 with the same starting point a.



 The idea is sketched in Fig. <ref> for two dimensions.


Take an irrational number, for example π, and when reaching a for the n-th time, follow the cycle T_i iff the n-th position in the binary expansion of π is i. The picture shows indeterministic behavior in two dimensions, which only happens with probability zero, it is however obvious that this construction can be lifted into three dimensions using Example <ref> as a gadget to obtain a similar chaotic path with positive probability. 

Fig. <ref> is also an example for a two-dimensional HGRN that leads to indeterminism independent by the starting point of the trajectory, and therefore to chaos, because every possible trajectory in 𝒩 will at some point pass a, meaning φ(a)=Φ(a)=𝒩.


This kind of chaos differs from the chaos that makes HGRNs Turing-complete. It depends on the choices made on the way rather than the restrictions of the network, and another trajectory starting in the same point might be periodic. To examine more precisely how the reachability problem and the identification of chaos relate to those, from now on we call a chaotic path that reaches a certain point an infinite number of times recurrent, and otherwise say that its chaos is proper.

The proof of Theorem <ref> relies on proper chaos, so reachability could still be decidable in recurrent chaos.

The following is a direct consequence of <cit.>:

A two-dimensional HGRN never admits proper chaos.


We now want to define a property that most HGRNs of natural origin share:


For a HGRN 𝒩 with two discrete levels for every gene and H=(V,E,ρ) a complete weighted directed graph with V={1,...,N} the genes of 𝒩 and ρ:E→ℚ, we say that 𝒩 is induced by G, iff for every α∈{0,1}^N we have that celerities are the weighted sums over the incoming edges

    c(α)_j=∑_i=1 
    α_i=1^Nρ(i,j)-∑_i=1 
    α_i=0^Nρ(i,j)

We call H the influence graph of 𝒩.


We want to pose the question whether our problems become easier when inducedness is given.

Note that induced HGRNs form a zero-set among all possible HGRNs.



It can be decided whether a HGRN is induced or not in P. Moreover, the influence graph is unique and can be computed in polynomial time.




When given all celerity vectors, the decisional as well as the computational problem are linear equation systems on rational domains and therefore solved efficiently. It remains to show that there can only be one solution.

Suppose there were two graphs H_1 and H_2 inducing the same HGRN, then for all dimensions i,j∈{1,...,n} we have that 


    ρ_1(i,j)=c([0]^i-1,1,[0]^n-i)_j-c(0,...,0)_j/2=ρ_2(i,j)

which gives us uniqueness.▪



The constellation in Example <ref> can be obtained by an induced HGRN with positive probability, for example by




    < g r a p h i c s >









§ CONCLUSION AND FURTHER QUESTIONS

We examined the reachability problem for HGRNs and showed in what way the size and dimension of the net influence chaos and indeterminism.

Further open questions are:



  1.) Are induced HGRNs still Turing-powerful?

  2.) Does Theorem <ref> still hold for induced HGRNs?

  3.) For H_𝒩 the set of all possible trajectories in 𝒩, how do the sets R_t:={h(t)| h∈ H_𝒩} behave for large t? With what probability and how fast do they approach a zero-set?

  4.) Can these sets R_t be computed efficiently?

  5.) Can HGRNs be expressed by induced HGRNs with the help of auxiliary dimensions?

Acknowledgment: We want to thank Klaus Meer for helpful discussion.




plain



