Boolean Networks with Intricate Updating Mode: Application to Genetic Regulatory Networks Involved in Familial Diseases

1. Introduction

Many familial diseases are due to accidents caused on the genome of the first pathological carrier in the family, by a problem of endogenous rewriting of the genetic message during meiosis, or by an exogenous influence (radiation, viral pathology, etc.). These accidents are generally nonsense point mutations that change only one amino acid or more important chromosomal recombinations of the genome, such as the translocations, causing often dramatic changes in the interactions between altered genes and the regulation RNAs (siRNAS, miRNAs, circRNAs, etc.) and proteins belonging to their epigenetic regulation graph. Figure 1 shows the graph of such interactions in the case of a regulation network centered on the ubiquitous transcription factor SP1, which plays a central role in many cellular pathways [1]. SP1 activates or represses the transcription of many other genes in response to physiological and pathological stimuli by binding with high affinity to GC-rich motifs of these genes involved in various processes such as cell apoptosis and growth as well as in differentiation and immune responses, cellular response to DNA damage, biological clock regulation, chromatin remodeling, oxidative stress response, etc. For this ubiquitous regulation role, SP1 is involved in many familial diseases such as Familial angioedema, Osteogenesis imperfecta and Biliary atresia. We will study in the following the attractors of the Boolean genetic networks regulating the main genes whose mutations cause these diseases. We will study the attraction basins of these attractors and show their robustness under five types of different perturbations (change of initial conditions, parameter values, interaction graph, transition function or updating clock) susceptible to provoke a loss stability of three types: trajectorial stability, asymptotic stability or structural stability. We will present the results in these two domains in Section 2, discuss these results and the chosen models in Section 3, describe the methods on both mathematical and numerical sides in Section 4 and conclude in Section 5 on the utility of robustness studies in the context of genetic regulation [2].

2. Results

In the following, we will examine in the framework of three applications the existence or absence of the stabilities associated to the last of the five types of perturbation described above, the change of updating mode, which is crucial in biology insofar as such an updating rule exists, because the exact functioning of the clock controlling the gene expression is unknown.

We will study the robustness of the genetic networks associated to three familial diseases [3,4,5,6]: the familial angioedema [7,8,9,10], the osteogenesis imperfecta [11,12,13] and the biliary atresia [14,15,16]

The genetic regulatory networks described in the following contain genes and interactions coming from the cited literature and from the reference databases [17,18,19].

2.1. Familial Angioedema

Figure 2 gives the network regulating the genes involved in angiœdema, a familial disease with a prevalence of about 1/50,000, with 2 types of expression of SERPING1, one with low (type I) and another with high (type II) expression level (the network has been provided initially by MetaCore^® [20] and checked after in the literature).

2.2. Osteogenesis Imperfecta Inducing Early Osteoarthritis

Figure 3 gives the network regulating the genes involved in osteogenesis imperfecta. All the interactions like Notch1 activating TGF through Pparga and M2 polarization come from [19,20,21,22].

2.3. Biliary Atresia

Biliary atresia regulation network comes from the top part of the network of Figure 4 centered on SP1 [14].

2.4. Dynamics of the Familial Angioedema Network

Recent data on the prevalence of familial angioedema come from [17]. They are in agreement with the expression data from Figure 5B and explain the choice of attractors 6 and 11 respectively for the type I and II of this hereditary disease. We can notice on Figure 5A that the physiologic attractor 3 is a fixed point (PF3) whose attraction basin is made of 93.43% of the possible states of the network. The pathologic attractors corresponding to the hereditary forms of the disease called angioedema are the fixed points PF6 and PF11, at which there is an absence of the gene SERPING1 in attractor 6 for type I and presence for type II. The relative size of their attraction basins relative size equals about 1‰ of the possible states of the network, and the observed prevalence of hereditary angioedema in the general population is reported to range from 1:10,000 to 1:150,000, without major sex or ethnic differences [23,24]. The ratio between the two relative sizes of pathologic attractors shows that type I represents 84% of the total as observed in real patients. E_attractor can be evaluated from Figure 5A where the number of genes considered is equal to 21 by supposing that circular RNAs are absent, and the attraction basin relative sizes are given on the line ABRS, then we have E_attractor ≈ 0.4520. Eventually, a previous work [25] gives an algebraic formula allowing for the calculation of the number of attractors of a Boolean network, where there are two tangential circuits (Figure 2C). In the top part of the network of Figure 2C, there are two such circuits, one positive of length five and one negative of length eleven. The attractor number predicted by the theory equals 11 and is the same as that obtained from the simulation of all trajectories of the network involved in angioedema from all possible initial conditions as summarized in Figure 5A.

2.5. Dynamics of the Osteogenesis Imperfecta Network

We propose to give precise results concerning the dynamics of a simplified version of the core of the network associated with the familial osteogenesis imperfecta (see Figure 3B and Figure 6), we will call network J in the following.

Figure 6 shows logical version of the transition function of the interaction graph involved in the osteogenesis imperfecta. By simulating in a parallel updating mode the iterations of the Boolean network having this transition motion from the initial condition where all the genes are inactive, we get from any initial configuration the following limit cycle of length 6 LC: (000000, 010000, 011000, 111000, 101000, 100000). The Table 1 exhibits all the possible dynamical behaviors of the network J depending on various initial configurations, showing the convergence to two attractors, the limit-cycles LC1 and LC2, the first of length 6 with a large basin of 56 configurations and the second of length 2 with a small basin of 2 configurations corresponding to pathologic states. Figure 7 and Table 2 corresponding to intricate updating show a unique limit-cycle LC of length 6 representing only normal states.

2.6. Dynamics of the Biliary Atresia Network

The number of attractors (7 attractors on Table 3) can be predicted by algebraic formula [25] allowing for the calculation of the number of attractors of a Boolean network with tangential circuits. For example, in case of 2 negative circuits of length 3 and 4, this number is 1. In the network of Figure 4, there are two such circuits, one negative of length five and one negative of length eleven. The attractor number predicted by the theory is 7 (see Table 4), the same as that obtained from the simulation of all trajectories of the network involved in biliary atresia from all possible initial conditions as summarized in Table 3.

The attractor limit-cycle LC2 of Table 1 or Table 3A corresponds to a high level of HRS (histone related sequences) and a low level of the expression of EGFR (inhibited by HRS on Figure 4B). Hence, it could correspond to the pathologic (BA) attractor of Table 3B, which shows an expression high for HSR and a low for EGFR (Figure 4B).

3. Discussion

The network J is related to familial angioedema was robust because, regardless of the block-sequential update mode used, fixed points will never appear, which is coherent with the empiric observation that most of genes involved the disease have a periodic expression due in part to the chromatin clock and circadian rhythms [1]. Some genes are updated systematically together under the control of the In order to achieve a systematic plan of numerical experience, we shall compare the block-sequential update modes of the Figure 5 to the intricate block sequential update mode, for which we consider that genes 3 and 4 are always active expressing their protein (Figure 6).

4. Methods

The definitions and theoretical results of this Section will follow previous works [26,27,28,29,30,31,32,33,34,35,36].

4.1. Kauffman Boolean Networks

We denote the whole configuration of a Kauffman Boolean network N of size n at time t by: x(t)={x_i(t)}_i=1,n, x(t) ∈ Ω={0,1}ⁿ, where Ω is the set of all possible Boolean configurations on N. Each node i of N represents a gene and is connected to k(i) other genes i₁,...,i_k(i) whose state is updated according to a specific order and a specific Boolean rule:

$$x_i\left(t+1\right)=F_i(x_{i_1},\dots ,x_{i_{k\left(i\right)}})$$

where F_i is a Boolean transition function and $x_{i_j}$represents the state of the gene i_j connected to i, where i_j may by equal to i itself. The Boolean function F_i is represented by a truth table which gives the output for each of the possible set of input values. For a function F_i with k(i) variables, there are 2^k(i) possible sets of input values, each yielding an output in {0,1}.

4.2. Interaction Graph

The interaction graph G(N) of the network N is represented by an oriented graph in which the n nodes correspond to genes and directed edges to the interactions between these genes. Each interaction is characterized by a logic operator w_ij (considered as an influence weight) which represents the influence that a source node of the interaction graph G(N), i.e., a gene j, exerts on an end node, i.e., a gene i, through the influence of its messenger RNA or of the protein the gene j expresses. More precisely, the w_ij is the identity represented by the sign + or the color green on the edge between j and I, or the negation represented by the sign – or the color red, depending on the fact that the gene j tends to respectively activate or inhibit gene i. There is no edge between j and i, if j has no influence on i.

4.3. Attractors

Under the influence of the transition functions F_i, the configuration x(t) of the network can evolve in Ω toward two possible asymptotic behaviors called attractors: a fixed configuration (also called fixed point) or a cycle of configuration (also called limit cycle). One of the challenges in a genetic regulation network is to determine the nature of its attractors (fixed point or limit cycle), which can have a great influence on its physiological or pathological character. For each attractor, we call attraction basin the set of all initial configurations of the network which lead after a finite number of iterations to this attractor. In the applications, we determine the possible attractors and their attraction basins (at least their size).

4.4. Frustration, Energy and Entropy

During the dynamics of the network, the states of two interacting genes j and i can change. If x_i(t) = x_j(t) when the weight between j and i is the identity, or if x_i(t) =ùx_j(t) when the weight between j and i is the negation, then the local frustration of the interaction between j and i is equal to 0. Else, this frustration equals 1. The global frustration of the network F(t) (or simply the network frustration) of the configuration x(t) is equal to the sum of its local frustrations. The network frustration is equal to the number of edges whose source x_j(t) and end x_i(t) states are in contradiction with their logic operator w_ij.

The kinetic energy E(t) of the network N at time t is defined by the following formula [2]:

$$E\left(t\right)={\left({\sum }_{i=1}^n{{(x}_i\left(t\right)-x_i(t-1))}^2\right)}^{{\displaystyle \frac{1}{2}}}/2$$

In the following, we will denote E(t) and call network energy the quantity E(t) = 4 E(t)². Network energy is equal to the number of nodes whose state changes between t and t+1.

The attractor entropy E_attractor is a measure of heterogeneity of the attractor landscape of the network dynamics defined by:

$$E_{attractor}=-{\sum }_{k=1,m}ABRS\left(A_k\right){log}_2\left[ABRS\left(A_k\right)\right],$$

(3)

where ABRS(A_k)=Card (A_k)/2ⁿ is the attraction basin relative size of the attractor A_k among the m attractors of the network.

4.5. Updating modes

An updating mode of the network is a rule describing in which order the state of the nodes has to be updated. There is three types of updating modes:

-

the block-sequential mode

It consists in choosing a partition of N in m disjoint subsets of nodes S₁,…,S_m with ${\cup }_{k=1}^m{S}_k=N$, which are updated sequentially, the nodes of each subset being updated parallelly. If each subset is a singleton, the updating mode is called sequential, the choice of the order being possibly random.

-

the block parallel mode

-

It consists in choosing a partition of N in m disjoint subsets of nodes S₁,…,S_m with ${\cup }_{k=1}^m{S}_k=N$_, which are updated parallelly, the nodes of each subset being updated sequentially. If m=1, S₁ = N, the updating mode is called parallel.

-

the block-intricate sequential of parallel modes

It consists in m non obligatory disjoint subsets of nodes S₁,…,S_m with ${\cup }_{k=1}^m{S}_k=N$, which are updated sequentially (respectively parallelly), the nodes of each subset being updated parallelly (respectively sequentially). An updating robustness study of network dynamics corresponds to consider all the possible updating modes and show which changes of states can occur when updating mode changes. The network can be robust for five types of perturbations (change of initial conditions, parameter values, interaction graph, transition function or updating clock) and three types of stability:

-

trajectorial (or Lyapunov) stability, which corresponds to the existence of a distance threshold respected between the ancient trajectory and the new after perturbation

-

asymptotic stability, which corresponds to the conservation of the number and nature of the attractors, even if the transient part of trajectories changes

-

structural stability, which corresponds to the conservation of the attraction basins in response to structural perturbations (change of interaction graph, transition function or updating clock).

4.6. Examples of Simulations

Simulations of Section 5.6. (Figure 8 and Figure 9) have been achieved using Boolnet [37].

In the parallel updating mode, the network K has 3 fixed points: 0 1 0 1 1 1 1 0 0, 1 0 1 0 0 0 0 1 1, 0 1 0 0 0 0 1 0 0 and two limit-cycles of length 3.

In the parallel updating mode, the network J has two attractors, the following limit-cycle of length 6: 000000 → 010000 → 011000 → 111000 → 101000 → 100000 → 000000 and the following attractor of length 2: 110000 → 001000.

5. Conclusions

We have defined a new update mode, the entangled block-sequential mode, in which the nodes of the Boolean automata network are updated differently, some being updated at each iteration and others being updated at lower frequencies (every other time, in the chosen example). In the application to Boolean networks of genetic regulation, in the absence of precise information on the mode of updating the states of the nodes of the network (i.e., the expression of the genes involved), it is necessary to examine the consequences of a change of mode on the dynamics of the network (possible change in the number, nature and size of the basins of attraction of its attractors). If the characteristics of the attractors are invariant in a change of update mode, we say that the network is robust. Consequently, the choice of a dynamic is less crucial for the practical consequences of the discrete modeling of the expression of the genes studied, in particular to explain the genesis of a pathology of familial origin involving a set of well-identified interacting genes, which is the case in the diseases taken examples, namely familial angioedema, osteogenesis imperfecta and biliary atresia.