Biological Networks Entropies: Examples in Neural Memory Networks, Genetic Regulation Networks and Social Epidemic Networks

Abstract

Networks used in biological applications at different scales (molecule, cell and population) are of different types: neuronal, genetic, and social, but they share the same dynamical concepts, in their continuous differential versions (e.g., non-linear Wilson-Cowan system) as well as in their discrete Boolean versions (e.g., non-linear Hopfield system); in both cases, the notion of interaction graph G(J) associated to its Jacobian matrix J, and also the concepts of frustrated nodes, positive or negative circuits of G(J), kinetic energy, entropy, attractors, structural stability, etc., are relevant and useful for studying the dynamics and the robustness of these systems. We will give some general results available for both continuous and discrete biological networks, and then study some specific applications of three new notions of entropy: (i) attractor entropy, (ii) isochronal entropy and (iii) entropy centrality; in three domains: a neural network involved in the memory evocation, a genetic network responsible of the iron control and a social network accounting for the obesity spread in high school environment.

1. Introduction

Biological networks are currently widely used to represent the mechanisms of reciprocal interaction between different biological entities and at different scales [1,2,3]. For example, if the entities are molecular, the networks can concern the genes and the control of their expression, or the enzymes and the regulation of the metabolism in living systems that possess them. At the other end of the life spectrum, at the population level, biological networks are used to model interactions between individuals, for example during the spread of contagious diseases, whether social or infectious [4,5]. Between these two molecular and population levels, the cellular level is also concerned by biological networks: for example, in the same neural tissue, cells such as neurons work together to give rise to emergent properties, such as those of the central nervous system related to memory [6] or motion control [7].

Among the tools used for quantifying the complexity of a biological network, the network entropies are ones of the most easy to calculate. They can be related to the network dynamics, i.e., calculated from the Kolmogorov-Sinaï formula when the dynamics is Markovian, or they can represent the richness of the network attractors or confinors [8], based on the relative weight of their attraction or confinement basins in the network state space. A final way to calculate the entropy of a random network is to use the invariant measure of its dynamics, giving birth to probability measures related to the characteristics of its nodes like centrality, connectivity, strong connected components size, etc. The literature on network entropies is abundant and from reference to recent papers it is possible to find definitions and calculation algorithms in [9,10,11,12,13,14,15,16,17,18,19,20,21,22].

Dynamical systems, discrete or continuous, share common mathematical notions like initial conditions and attraction basin, final conditions and attractor, Jacobian graph, stability and robustness. Moreover, in an energetic view, there are in both cases, two extreme systems, i.e., on one hand, conservative systems with an energy constant along trajectories (e.g., Hamiltonian systems) and, on the other hand, dissipative systems with an energy decreasing along trajectories (e.g., potential systems). For both Hamiltonian and potential systems, it is possible to define at each point x of their state space E, a Jacobian matrix J(x) made of the partial derivatives of the flow of the network dynamics, whose associated graph having J(x) as incidence matrix, is called the interaction graph G(x). The description of the dynamics in terms of energy conservation or dissipation is identical whether the system is time discrete or continuous.

Common definitions of the concepts of attraction basin and attractor will be recalled in Section 2 Material and Methods, in which we give also some examples of Hamiltonian and potential systems by making explicit the energy functions underlying their dynamics, allowing for interpreting these functions in the precise context of different applications. For example, the energy functions can be related to the notions of frustration in statistical mechanics [23] or of kinetic energy, which can be calculated in discrete systems by using an analogue of the continuous velocity [24]. In Section 3, we present the Results through applications in three domains, neural (memory evocation based on the synchronization of neural networks), genetic (regulation of the iron metabolism) and social (dynamics of the obesity spread in a college environment). We define the thermodynamical notions of attractor entropy, dynamic (or evolutionary) entropy and entropy centrality for interpreting the results, and we relate the dynamic entropy to the robustness (or structural stability) of the dynamical systems presented. After, in Section 4, the Discussion presents the interest of a thermodynamic approach in social networks in order to simulate and evaluate the efficiency of public health policies, and eventually, the Section 5 gives the Conclusion and Perspectives of the paper.

2. Material and Methods

2.1. The Dynamical Tools

Several definitions of attractor and stability have been proposed since 50 years. Here we use that of Cosnard and Demongeot [25,26], which covers both continuous and discrete cases. More specific approaches of the notion of attractors especially for network dynamics can be found in [3,27,28,29].

2.1.1. Definition of the Notion of Attractor

Let consider a trajectory ${\left\{f(a,t)\right\}}_{t\in \mathrm{T}}$ starting at the point a of the state space $E\subset {\mathbb{R}}^n$ provided with a dynamic flow f and a distance d. $L\left(a\right)$, called the limit set of the trajectory starting in a, is the set of accumulation points of the trajectory ${\left\{f(a,t)\right\}}_{t\in \mathrm{T}}$ as time t goes to infinity in the time set T (either discrete or continuous) such as:

$$L\left(a\right)=\left\{y\in E;\forall \epsilon >0,\forall t\in \mathrm{T},\exists s(\epsilon ,t)>t/d\left(f\right(a,s(\epsilon ,t)),y)<\epsilon \right\}$$

where $s(\epsilon ,t)$ is one time greater than t for which the trajectory passes into a ball of radius $\epsilon$ centred in y.

If the initial state a lies in a set A, $L\left(A\right)$ is the union of all limit sets $L\left(a\right)$, for a belonging to A:

$$L\left(A\right)={\cup }_{a\in A}L\left(a\right)$$

Conversely, $B\left(A\right)$ is the set of all initial conditions a, whose limit set $L\left(a\right)$ is included in A (Figure 1). $B\left(A\right)=\left\{y\in E;L\left(y\right)\subset A\right\}$ is called the attraction basin of A.

(i)

$A=LoB\left(A\right)$, where the composed operator $LoB$ is obtained by applying the basin operator B and then the limit operator L to the set A

(ii)

There is no set A’ containing strictly A, shadow-connected to A and verifying (i)

(iii)

There is no set A” strictly contained in A verifying (i) and (ii)

The definition of the “shadow-connectivity” connecting in E the subsets $F_1$ and $F_2$, lies on the fact that there exists a “shadow-trajectory” between $F_1$ and $F_2$ (Figure 1 Top right, where $F_1=\left\{x\right\}$ and $F_2=\left\{y\right\}$). The notion of “shadow-trajectory” has been defined first in [30]:

x in $F_1$ and y in $F_2$ belong to the same shadow-trajectory between $F_1$ and $F_2$ if, for any $\epsilon >0$, there is an integer $n\left(\epsilon \right)$ and $n\left(\epsilon \right)$ times $t(1,\epsilon ),\dots ,t\left(n\right(\epsilon ),\epsilon )$ and states $y(1,\epsilon ),\dots ,y\left(n\right(\epsilon ),\epsilon )$, such as all Hausdorff distances between the successive parts of the shadow-trajectory are less than $\epsilon$:

$$d\left(f\right(x,t(1,\epsilon )),y(1,\epsilon \left)\right)\le \epsilon ,\dots ,d\left(f\right(y(k-1,\epsilon ),t(k-1,\epsilon )),y(k,\epsilon \left)\right)\le \epsilon ,\dots ,d\left(f\right(y\left(n\right(\epsilon ),\epsilon ),t\left(n\right(\epsilon ),\epsilon )),y)\le \epsilon$$

The above definition of attractor is available for all known dynamical systems, in agreement with the common sense and is meaningful for both discrete or continuous time sets T, if the dynamical system is autonomous with respect to time, contrary to other definitions proposed since 1970 such as in [31,32,33,34], all less general.

2.1.2. Potential Dynamics

A dynamical system is potential (or gradient) if its velocity is equal to the gradient of a scalar potential P defined on the state space E: $\partial x\left(t\right)/\partial t=-gradP=-\partial P/\partial x$ where gradient vector equals $\partial P/\partial x=(\partial P/\partial x_1,\dots ,\partial P/\partial x_n)$ and state at time t is $f(x\left(0\right),t)=x\left(t\right)=(x_1\left(t\right),\dots ,x_n\left(t\right))$ [35]. The system is dissipative, because P decreases along trajectories until its attractors located on its minima.

2.1.3. Hamiltonian Dynamics

A dynamical system is Hamiltonian if its velocity is tangent to the projection on E of the contour lines of the surface representative of an energy function H on E: $\partial x\left(t\right)/\partial t=tangH\left(x\right(t\left)\right)$, where $tangH\left(x\right(t\left)\right)$ is the vector tangent at $x\left(t\right)$ to the projection on E of the equi-altitude set of altitude $H\left(x\right(t\left)\right)$. If the dimension of the system is two, then $tangH$ is equal to $(\partial H/\partial x_2,-\partial H/\partial x_1)$ [36]. The system is said conservative, because the energy function H is constant along trajectories. An example of conservative system will be given in Section 4.2 (Equation (16)).

2.1.4. Potential-Hamiltonian Dynamics

A dynamical system is potential-Hamiltonian if its velocity can be decomposed into two parts, one potential and one Hamiltonian: $\partial x\left(t\right)/\partial t=-gradP+tangH$ [37]. If the set of minima of P is a contour line of the surface H on E, then its strong shadow-connected components (i.e., the sets in which any pair of states is shadow-connected in both senses) are the attractors of the system.

2.2. Examples of Dissipative Energy

In a gradient system, the potential energy is dissipative, because it decreases along the trajectories and, more, its gradient is exactly equal to the system velocity. The existence of a morphogenetic dissipative energy has been speculated by the Waddington’s school of embryology [38] as an ontology rule governing the development of living systems and we give in Section 2.2.2 an example of such a dissipative energy.

2.2.1. Discretization of the Continuous Potential System with Block-Parallel Updating

Let us consider a discrete Hopfield dynamics where $x_i\left(t\right)=1$, if the gene i is expressing its protein, else $x_i\left(t\right)=0$. The probability that $x_i\left(t\right)=1$, knowing the states in an interaction neighbourhood $V_i$ of i is given by:

$$P\left(\left\{x_i\left(t\right)=1\right\}\mid \left\{x_j(t-1),j\in V_i\right\}\right)=exp\left(\left(\sum _{j\in V_i}{w}_{ij}{x}_i(t-1)-\theta \right)/\mathrm{T}\right)$$

(1)

where $w_{ij}$ is a parameter which represents the influence (inhibitory, if $w_{ij}<0$; activatory, if $w_{ij}>0$ and null, if $w_{ij}=0$) the gene j from the interaction neighbourhood $V_i$ is exerting on the gene i, and $\theta$ is a threshold corresponding to the external field in statistical mechanics and here to a minimal value of activation for the interaction potential ${\sum }_{j\in V_i}{w}_{ij}{x}_i$. T is a temperature allowing for introducing a certain degree of indeterminism in the rule of the change of state. In particular, if $T=0$, the rule is practically deterministic and, if H denotes the Heaviside function, we have:

$$\begin{array}{c}\hfill x_i\left(t\right)=1, \mathrm{if}\sum _{j\in V_i}{w}_{ij}{x}_i(t-1)>\theta \\ \hfill x_i\left(t\right)=0, \mathrm{if}\sum _{j\in V_i}{w}_{ij}{x}_i(t-1)<\theta \end{array}\}{x}_i\left(t\right)=H(\sum _{j\in V_i}{w}_{ij}{x}_i(t-1)-\theta )$$

(2)

$$P(\left\{x_i\left(t\right)=1\right\}\mid \left\{x_j(t-1),j\in V_i\right\})=1/2,\mathrm{if}\sum _{j\in V_i}{w}_{ij}{x}_i(t-1)=\theta$$

Once this rule defined, three graphs and their associated incidence matrices are important for making more precise the network dynamics. These incidence matrices are:

• interaction matrix W: $w_{ij}$ represents the action of the gene j on gene i. W plays the role of a discrete Jacobian matrix. We can consider the associated signed matrix S, with $s_{ij}=1$, if $w_{ij}>0$, $s_{ij}=-1$, if $w_{ij}<0$ and $s_{ij}=0$, if $w_{ij}=0$,
• updating matrix U: $u_{ij}=1$, if j is updated before or with i, else $u_{ij}=0$,
• trajectory matrix F: $f_{bc}=1$, where b and c are two states of E, if and only if $b=f(c,1)$, else $f_{bc}=0$.

These matrices can be constant or can depend on time t:

• in the case of W, the most frequent dependence is called the Hebbian dynamics: if the vectors ${\left\{x_i\left(s\right)\right\}}_{s<t}$ and ${\left\{x_j\left(s\right)\right\}}_{s<t}$ have a correlation coefficient ${\rho }_{ij}\left(t\right)$, then $w_{ij}(t+1)=w_{ij}\left(t\right)+h{\rho }_{ij}\left(t\right)$, with $h>0$, corresponding to a reinforcement of the absolute value of the interactions $w_{ij}$ having succeeded in inhibiting or activating their target gene i: in case where, for $s<t$, $x_j\left(s\right)$ remained equal to one, that leads to increase the $x_i\left(s\right)$’s, if the $w_{ij}\left(s\right)$’s were positive, and conversely to decrease the $x_i\left(s\right)$’s , if the $w_{ij}\left(s\right)$’s were negative,
• in the case of U, we can have an autonomous (in time) clock based on the behaviour of r chromatin clock genes having indices $1,\dots ,r$, with three possible behaviours:
  • if $y\left(t\right)={\prod }_{i=1,\dots ,r}{x}_i\left(t\right)=1$, then the rule (1) is available
  • if $y\left(t\right)=0$ and ${\sum }_{s=t,\dots ,t-c}y\left(s\right)>0,x(t+1)=x\left(s^{*}\right)$, where $s^{*}$ is the last time between $t-c$ and t, where $y\left(s^{*}\right)=1$
  • if $y\left(t\right)=0$ and ${\sum }_{s=t,\dots ,t-c}y\left(s\right)=0$, then $x(t+1)=0$ (by exhaustion of the pool of expressed genes).

This dynamical system ruling the evolution of U remains autonomous, but the updating is state dependent. If all $u_{ij}$ equal one, the updating is called parallel; if there is a sequence $i_1,\dots ,i_n$ such as $u_{kk+1}=1$ for any k such as $1\le k\le n-1$, the other $u_{ij}$ being equal to zero, then the updating is called sequential. There exist cases other than these extreme cases called block-sequential (the updating is parallel in a block and the blocks are updated sequentially). A case of updating which is more realistic in genetics is called block parallel: there are blocks made of genes sequentially updated, and these blocks are updated in parallel. Interactions between these blocks can exist (i.e., there are $w_{ij}\ne 0$, with i and j belonging to two different blocks), and because the block sizes are different, the time at which the stationary state is reached in a block is not necessarily synchronized with the corresponding times in other blocks: the synchronization depends on the intra-block as well as on the inter-block interactions. All the genes of a block can reach their stationary state long before the genes of another block. This behaviour explains in particular that the state of the genes of the block corresponding for example to the chromatin clock is highly depending on the state of the genes in other blocks connected to them (and for example acting as transcription factors of the clock genes), and reciprocally.

An example of such a block-parallel dynamics is given by a short block made of the nodes 4 and 5 on Figure 2 Top middle of chromatin clock genes reaching rapidly its limit-cycle and entraining other blocks of metabolic (nuclear or mitochondrial) genes. The second block of genes is made of the morphogens responsible of the apex growth of a plant and of the auxin gene inhibiting the morphogens (node 1) [39,40]. When the first block reaches its limit-cycle (possibly entrained by external “Zeitgebers” related to an exogenous seasonal rhythm), it blocks the dynamics of the second block, liberating the first and after, the second axillary bud morphogens (nodes 2 and 3) only when the apex morphogens passed a sufficient time in state 1, corresponding to a sufficient growth of apex for limiting the diffusion of auxin, then cancelling its inhibitory power. We can easily iterate this dynamical behaviour for other successive blocks related to the inferior buds.

More precisely, we can implement the dynamics of the morphogenesis network in the following way; let us suppose that:

(i)

there are two chromatin clock genes involved in a regulon, i.e., the minimal network having only one negative circuit and one positive circuit reduced to an autocatalytic loop

(ii)

there are three morphogens corresponding to the apex and to axillary buds involved as nodes of a 3-switch network fully connected, with all interactions negative except three autocatalytic loops at each vertex

(iii)

there are three inhibitory interactions from the auto-catalysed node of the regulon on the nodes of the 3-switch, as indicated on Figure 2 Top left.

The Figure 2 Top shows the three graphs defining the morphogenesis network dynamics, i.e., the interaction (left), updating (middle) and trajectory (right) ones. The dynamics is supposed ruled by a deterministic Hopfield network, with temperature and threshold equal to zero and interaction weights $w_{ij}$ equal to one (if positive), −1 (if negative) or zero. An attractor limit-cycle of this dynamics is shown on Figure 2 Top right, showing two intricate rhythms, one of period four corresponding to the chromatin clock dynamics, embedded in a rhythm of period 12 corresponding to the successive expression of the apex and axillary buds morphogens. The Figure 2 Bottom shows other updating graphs, less realistic than the block-parallel for representing the action of the chromatin clock regulon.

2.2.2. Discrete Lyapunov Function (Neural Network)

A logic neural network has local transition rules of type ⊕, ∨ and ∧ (Figure 3). Applying these local rules ∧ (“and”), ∨ (inclusive “or”) et ⊕ (exclusive “or”), we can define on each configuration $x\left(t\right)$ of the logic neural network of the Figure 3, a Lyapunov function L from the global frustration of order four $D_4$, decreasing along trajectories and vanishing on the attractors of the network dynamics:

$$L\left(x\left(t\right)\right)=\left[D_4\left(x(t+1)\right)+D_4(x\left(t\right)\right]/64$$

with $D_4\left(x\left(t\right)\right)={\sum }_{i=1,\dots ,n}{D}_{4,i}\left(x\left(t\right)\right)$, where $D_{4,i}$ is the local frustration of order four, also equal to the discrete kinetic energy $E_{4,i}\left(x\left(t\right)\right)$ between $x_i\left(t\right)$ et $x_i(t+4)$ multiplied by 32:

$$D_{4,i}\left(x\left(t\right)\right)/32={(x_i(t+4)-x_i\left(t\right))}^2/2{(t+4-t)}^2=E_{4,i}\left(x\left(t\right)\right)$$

In this example of a discrete neural network ruled by logic Boolean functions applied for each node to its entries (proposed in [41,42] as a genetic network), the dynamical behaviour is guided by an energy function deriving from the discrete kinetic energy, built with the discrete equivalent of the velocity [24]. We can indeed write the Lyapunov function of the trajectory $L\left(x\right(t\left)\right)$ as the mean global discrete kinetic energy calculated between times t and $t+5$:

$$L\left(x\left(t\right)\right)=\left[D_4\left(x(t+1)\right)+D_4(x\left(t\right)\right]/64=[E_4\left(x(t+1)\right)+E_4\left(x\left(t\right)\right)]/2$$

3. Results

3.1. Attractor Entropy and Neural Networks

The attractor entropy $E_{attractor}$ measures the heterogeneity of the attractor landscape on the state space E of a dynamical system. It can be evaluated by the quantity:

$$E_{attractor}=-\sum _{k=1,\dots ,m}RABS\left(A_k\right)Log\left[RABS\left(A_k\right)\right]$$

(3)

where $RABS\left(A_k\right)$ is equal to the attraction basin size of the kth attractor $A_k$ among the m attractors of the network dynamics, divided by $2^n$, n being the number of genes of the network. We will give an example of use of this notion in the case of a continuous neural network, the Wilson-Cowan oscillator, which describes the biophysical interactions between two populations of neurons:

$$dx/dt=-x/{\tau }_x+tanh(\lambda x-\theta )\text{–}tanh(\lambda y-\theta ),dy/dt=-y/{\tau }_y+tanh(\lambda x-\theta )+tanh(\lambda y-\theta )$$

(4)

The variables x and y represent the global electrical activity of the two populations, respectively excitatory and inhibitory, whose interaction is described through a sigmoid function parametrized by $\lambda$, which controls its stiffness, accounting for the synaptic response to an input signal, which is non-linear when $\lambda$ is large (>>1), depending on the value of the input, greater or less than a potential barrier (here $\theta /\lambda$), and is an affine function of the input signal when $\lambda$ and $\theta$ are small (<<1). The parameters ${\tau }_x$ and ${\tau }_y$ refer to the membrane potential characteristics of the neuronal populations having respective activity x and y. If $\theta =0$ and ${\tau }_x={\tau }_y=\tau$, the Wilson-Cowan oscillator presents a Hopf bifurcation when $\lambda \tau$ crosses the value one [43], moving from a dynamics involving only an attracting fixed point to a dynamics comprising a repulsor at the origin of E, surrounded by an attractor limit-cycle C. If T denotes the period of the limit-cycle C, let us consider the point $h=f(O,h)$ of C reached after a time (or phase) h on the trajectory f starting on C at a reference point called the origin $O=f(O,0)$. In Figure 4 Top, $T=2\pi$ and this trajectory on C defines a one-to-one map F between the interval $[0,T]$ and C. We introduce now the notion of isochron $I_h$ of phase h in the attraction basin $B\left(C\right)$ as the set of points x of $B\left(C\right)$ such as $d\left(f\right(x,t),f(h,t\left)\right)$ tends to zero, when t tends to infinity. $I_h$ extends the map F to all the points of B(C). If we consider as new state space $I_h$, a new time space $T_h=\{t\in T;t=h+kT,k\in \mathbb{N}\}$ and the restriction of the flow f on $I_h$, such as: $\forall x\in I_h, f(x,t)\in I_h, h=li{m}_{t\in T_h;t\to \infty }f(x,t)$, we have: $B\left(C\right)={\cup }_{h\in C}{I}_h$.

Then, h can be considered as an attractor for f, whose basin is $I_h$. Let us decompose now the interval $[0,T]$ in n equidistant sub-intervals ${\left\{[kT/n,(k+1)T/n]\right\}}_{k=0,\dots ,n-1}$. If we consider now the time space $[kT/n,(k+1)T/n]$, the subset of C denoted $A_k=F\left([kT/n,(k+1)T/n]\right)$ is an attractor with an attraction basin equal to $B_k={\cup }_{h\in \left(\right[kT/n,(k+1)T/n\left]\right)}{I}_h$, for the flow restrained to $B_k$, and we have: $B\left(C\right)={\cup }_{k=0,\dots ,n-1}{B}_k$. Eventually, we can define the entropy related to these attractors in E, called isochronal entropy of degree n, which is just the attractor entropy calculated for the dynamics having as flow $f_h$:

$$E_{attractor}^n=-\sum _{k=1,\dots ,n}RABS\left(A_k\right)Log\left[RABS\left(A_k\right)\right]$$

(5)

If all the isochron basins have same size, then $E_{attractor}^n=Logn$. If not, the quantity $E_{attractor}^n$ reflects a spatial heterogeneity of the dynamics: certain $B_k$ can be wide (where the flow f is rapid inside) and other thin (where the flow f is slow inside). In the example of the Wilson-Cowan oscillator, when ${\tau }_x={\tau }_y=1$, $\theta =0$ and $\lambda$ is growing from zero to infinity, the isochronal entropy is diminishing from $Logn$ to $Log4$.

It is possible to change the value of the period T of the Wilson-Cowan oscillator by tuning the parameters ${\tau }_x$ and ${\tau }_y$ (Figure 4 Bottom right). Consider now values of the parameters, ${\tau }_x={\tau }_y=\tau =1$ and $\lambda =1.1$, close to the bifurcation in the parameter space $\left(\lambda \tau \right)$ and superimpose the isochron landscape on a map in false color giving the asymptotic phase shift $\left|\right|h^{\prime }-h\left|\right|$ between a trajectory starting at a point x on the isochron $I_h$ and a trajectory starting, after an instantaneous perturbation e, at the point x+e belonging to the isochron $I_{h^{\prime }}$ (Figure 4 Bottom left). The maximum phase shift of the Wilson-Cowan oscillator is inversely proportional to the intensity of this instantaneous perturbation, which means that a population of oscillators can then be synchronized by a sufficiently intense perturbation. Indeed, isochrons are spirals, which diverge from one to another by moving away from the attractor limit-cycle. The phenomenon of memory evocation after a sensory stimulation relies on the synchronizability of a population of identical (but dephased) oscillators whose successive states of their limit cycle store an information [44,45]. This evocation is the more effective since all oscillators have the same phase after stimulation, which is the case if the corresponding perturbation done at different phases ${\phi }_1$ and ${\phi }_2$ on the limit-cycle leads the oscillator between two successive isochrons (Figure 4 Bottom left). By using either the discrete version of the Wilson-Cowan system, namely the Kauffman-Thomas cellular automaton (the deterministic version of the Hopfield model), or its continuous version (Figure 4), it is easy to see that the potential part of the system is responsible for the spacing between isochrons away from the limit-cycle, causing the decrease of the isochronal entropy, and that the Hamiltonian part is responsible for their spiralization (Figure 4 Top).

By playing with the potential-Hamiltonian distribution [46], it is therefore possible, if the system is instantly perturbed, to increase the efficiency of the post-perturbation synchronization (due to the potential return to the limit-cycle) or to allow for desynchronizing (due to the Hamiltonian spiralization), which is necessary to avoid the perseveration in a synchronized state (observed in many neuro-degenerative pathologies).

3.2. Dynamic Entropy and Genetic Networks

We define first the functions energy U, frustration F and dynamic entropy E of a genetic network N with n genes in interaction [47,48,49,50,51,52,53,54,55,56,57,58,59].

$$\forall x\in \mathsf{\Omega },U\left(x\right)=\sum _{i,j\in \left\{1,n\right\}}{\alpha }_{ij}{x}_i{x}_j=Q_{+}\left(N\right)-F\left(x\right),$$

(6)

where x is a configuration of gene expression ($x_i=1$, if the gene i is expressed and $x_i=0$, if not), $\mathsf{\Omega }$ denotes the set of all configurations of gene expression (i.e., the hypercube ${\left\{0,1\right\}}^n$) and ${\alpha }_{ij}$ is the sign of the interaction weight quantifying the influence the gene j exerts on the gene i: ${\alpha }_{ij}=+1$ (resp. −1), if j is an activator (resp. inhibitor) of the expression of i, and ${\alpha }_{ij}=0$, if there is no influence of j on i. $Q_{+}\left(N\right)$ is the number of positive edges of the interaction graph G of the network N, whose incidence matrix is given by $A={\left({\alpha }_{ij}\right)}_{i,j\in \left\{1,n\right\}}$. $F\left(x\right)$ denotes the global frustration of x, i.e., the number of pairs of genes $(i,j)$ for which the values of $x_i$ and $x_j$ are contradictory with ${\alpha }_{ij}$, the sign of the interaction of j on i:

$$F\left(x\right)=\sum _{i,j\in \left\{1,n\right\}}{F}_{ij}\left(x\right),$$

where $F_{ij}$ is the local frustration of the pair $(i,j)$ equal to:

$$\left\{\begin{array}{c}1,\mathrm{if}{\alpha }_{ij}=1,\mathrm{with}{x}_j=1-x_i,\mathrm{or}{\alpha }_{ij}=-1,\mathrm{with}{x}_j=x_i\hfill \\ 0,\mathrm{elsewhere}\hfill \end{array}\right.$$

The dynamics of the network can be defined by a deterministic transition operator like the following threshold Boolean automaton with a majority rule, we will call in the following deterministic Hopfield transition:

$$x_i\left(t\right)=H(\sum _{j\in N_i}{w}_{ij}{x}_j\left(t\right)-{\theta }_i),$$

(7)

where H is the Heaviside step function defined by: $H\left(y\right)=1$, if $y>0$; $H\left(y\right)=0$, if $y\le 0$, $w_{ij}$ is the interaction weight of j on i, ${\theta }_i$ an activation threshold (whose exponential is called tolerance h: $\theta =logh$) and $N_i$ is the set of genes j such as $w_{ij}\ne 0$. The

Table 1. Recapitulation of the attractors of the iron regulatory network of Figure 5, for the parallel updating mode, with the list of expressed (in state 1) and not expressed (in state 0) genes and their relative attraction basin sizes, obtained by simulating in parallel updating mode a threshold Boolean network having a constant absolute value of non zero weights equal to one and a threshold equal to zero, from all the possible initial configurations.

Order	Gene	Fixed Point	Fixed Point 2	Limit-Cycle 1			Limit-Cycle 2
1	TfR1	0	0	0	0	1	1	0
2	FPN1a	0	0	0	0	0	0	0
3	C-Myc	0	1	0	0	0	0	0
4	Notch	0	0	1	1	1	1	1
5	GATA-3	0	0	1	1	1	1	1
6	IRP	0	0	0	1	0	0	1
7	Ft	0	0	0	0	0	0	0
8	Fe	0	0	0	0	1	0	1
9	MiR-485	0	0	0	0	0	0	0
10	CiRs-7 anti-sense	0	0	0	0	0	0	0
	Relative Attraction Basin Size	512/1024	256/1024	216/1024			40/1024

shows the parallel simulation of such a Boolean automaton, with $|w_{ij}|=1$ and ${\theta }_i=0$. This deterministic dynamics is just a particular case of the random dynamics defined by the transition analogue to the Equation (7):

$$m_{ij}=P\left(\left\{x_i\left(t\right)=1\right\}|\left\{x(t-1)=y\right\}\right)=exp\left[\sum _{j\in \{1,n\}}{w}_{ij}{y}_j/T\right]=P_{i,1}^{\left\{i/y_i=1\right\}},$$

(8)

when threshold and temperature equal to zero.

$m_{ij}$ is the probability of expression of the gene i at time t, knowing the configuration at time $t-1$ is equal to y and from the ${\left\{m_{ij}\right\}}_{i,j\in \{1,n\}}$ we can easily calculate the Markov transition matrix of the dynamics (8) supposed to be sequential [48,60], i.e., with a sequential updating of the genes in a predefined order.

Let $M=\left(M_{xy}\right)$ denote the Markov matrix giving the transition probabilities, defined by the rule (8) and an updating mode, between the configurations x and y of $\mathsf{\Omega }$, and $\mu =\left({\mu }_x\right)=(\mu {\left(\left\{x\right\}\right)}_{x\in \mathsf{\Omega }}$ be its stationary distribution on $\mathsf{\Omega }$. The dynamic entropy E can be calculated as:

$$E=-\sum _{x,y\in \mathsf{\Omega }}{\mu }_x{M}_{xy}log{M}_{xy}.$$

In the sequential updating mode, where the updating order of the nodes is the integer ranking $1,\dots ,n$, we have, by denoting $I=\{1,\dots ,i-1\},N\backslash I=\left\{i,\dots ,n\right\}$ and identifying x with the set of the indices i such that $x_i=1$:

$$M_{xy}=\prod _{i=1,\dots ,n}\left[P_{i,1}^{\left[x\cap (N\backslash I)\right]\cup \left[y\cap I\right]}{I}_{\left\{i\in y\right\}}+P_{i,0}^{\left[x\cap (N\backslash I)\right]\cup \left[y\cap I\right]}{I}_{\left\{i\notin y\right\}}\right]$$

and $\mu$ is classically the Gibbs measure [48,60] defined by:

$$\forall x\in \mathsf{\Omega },{\mu }_x=exp((\sum _{i\in x,j\in N_i}{w}_{ij}{x}_i{x}_j-{\theta }_i)/T)/Z,$$

where $Z={\sum }_{y\in \mathsf{\Omega }}exp(({\sum }_{j\in y,k\in N_j}{w}_{jk}{y}_j{y}_k-{\theta }_j)/T)$. When $T=0$, $\mu$ is concentrated on the $m(\le 2^n)$ attractors of the deterministic dynamics, then, its entropy is defined by:

$$E_{\mu }=E_{attractor}=-\sum _{k=1,\dots ,m}RABS\left(A_k\right)logRABS\left(A_k\right),$$

(9)

where m is the number of the final classes of the Markov transition matrix related to the dynamics (8). When $T=+\infty$, $\mu$ is scattered uniformly over $\mathsf{\Omega }$ and $E=nlog2$. An example of calculation of Relative Attraction Basin Sizes is given in the Table 1. We will estimate E between $T=0$ and $T=+\infty$ from the attractor entropy $E_{attractor}$ [51] by using the following approximate equality:

$$E\approx log{2}^n-E_{attractor}$$

(10)

E serves as robustness parameter, being related to the capacity the genetic network has to return to the equilibrium after endogenous or exogenous perturbations [53,59]. We can more generally prove that the robustness of the network, supposed to be sequentially updated, decreases when the variance of the frustration F of the network increases, because of the formula: $\partial E/\partial Log\left(\right|w\left|\right)=-VarF$, where $\left|w\right|$ is the absolute value supposed to be the same for all the non zero interaction weights between the genes of the network [51].

The regulatory network controlling the iron (Fe) metabolism is given in Figure 5. The central proteins controlling the iron metabolism are iron regulatory proteins IRP1 and IRP2 represented in Figure 5 by a unique entity called IRP, which can be active or inactive, active form only having a regulatory role. In absence of iron, IRP shows a mRNA-binding activity on a specific RNA motif called IRE, and when iron concentration is sufficiently high, it looses its mRNA-binding activity on IRE motif. The iron node (Fe) in Figure 5 represents the cellular iron available for various cellular processes and the node ferritin (Ft) represents a protein capable of internalize thousands of iron ions in a given Ft protein.

The complex mechanism of iron loading into, and release from, ferritin is taken into account here in a very simplified way: ferritin mRNA contains two IRE motifs in its 5’-UTR regions, and IRP consequently inhibits its translation. These IRE motifs are also inhibited by micro-RNAs like miR485, itself inhibited by the cyclic ARN ciRs7. IRP is also directly inhibited by ciRs7-antisense and activated by GATA-3, itself controlled by the positive circuit Notch/c-Myc The transferrin receptor (TfR1) located at the cellular surface allows iron uptake. TfR1 mRNA contains five IRE motifs in its 3’-UTR region, hence, the density of TfR1 receptors at the cell surface is correlated with the IRP activity. Eventually, ferroportin (FPN1a) is an iron exporter located at the cellular surface and contains one IRE motif in its 5’-UTR mRNA region, thus its translation is inhibited by IRP.

By studying the dynamics of the iron regulatory network (supposed to be Boolean, that is with only two states 1 and 0 for each entity, the state 1 corresponding to the proportion of its active form over a certain threshold), it is possible to calculate its asymptotic states (i.e., observed when the time goes to infinity) obtained for any initial condition in attraction basins (Table 1). From this dynamics study, we can then easily calculate $E_{attractor}=2.59412$, hence $E=4-2.59412=1.40588$, which corresponds to a relatively low robustness, the maximum being equal to four.

The notion of robustness is related to the capacity to diminish the consequences, on the entropy E, of endogenous or exogenous variations affecting the common absolute value $\left|w\right|$ of the non zero interaction weights of the genetic network. This notion of robustness is closed to that of structural stability in dynamical systems. If the variability of the global frustration F of the network is important, then its robustness diminishes. It is for example the case if a gene belonging to a positive circuit [61,62]. of a strong connected component of the interaction graph of the network is randomly inhibited. In this case, the inhibited gene behaves as a sink (i.e., without outing interaction) of state 0, which can enter in contradiction with that predicted by the states and interactions of its neighbors, increasing the variance of F and consequently diminishing the dynamic entropy E.

3.3. Centrality of Nodes and Social Networks

There are four classical types of centrality in a graph (Figure 7). The first is the betweenness centrality [63]. It is be defined for a node k in the interaction graph G as follows:

$$C_B\left(k\right)=\sum _{i\ne j\ne k\in G}\frac{{\beta }_{ij}\left(k\right)}{{\beta }_{ij}},$$

(11)

where ${\beta }_{ij}\left(k\right)$ is the total number of shortest paths from node i to node j that pass through node k, and ${\beta }_{ij}={\sum }_{i\ne j\in G}{\beta }_{ij}\left(k\right)$.

Let consider now a social network with a threshold Boolean dynamics following the Equation (7), containing overweight or obese (state 1) and normal (state 0) individuals as nodes. By concentrating the pedagogic efforts on the betweenness central nodes overweight or obese with tolerance 1, which are “hubs” of the social network, it results a rapid decrease of the number of overweight or obese, more important if the educative effort has been done on individuals with tolerance 1 than with tolerance 0 (Figure 6 Bottom right), as well as a decrease of their mean weight (Figure 6 Bottom left). The second corresponds to the degree centrality, which can be defined from the notions of out-, in- or total degree of a node i, corresponding to the number of arrows of the interaction graph, respectively outing, entering, or both, connected to i. For example, the in-degree centrality is defined by:

$$C_i^{indeg}=\sum _{j=1,\dots ,n;j\ne I}{a}_{ij}/(n-1),$$

(12)

where the general coefficient of the adjacency matrix A of the graph, denoted $a_{ij}$, equals one if there is a link between j and i, and else zero.

The third type of centrality is the closeness centrality. The closeness is the reciprocal of the average fairness, which is defined by averaging over all nodes j of the network else than i a distance between i and j, equal to ${\sum }_{j=1,\dots ,n;j\ne i}d(i,j)/(n-1)$, then:

$$C_i^{clo}=(n-1)/[\sum _{j=1,\dots ,n;j\ne i}d(i,j)]$$

(13)

where the chosen distance is the length of the closest path between i and j. The last type of classical centrality is the spectral centrality or eigen-centrality, which takes into account that the neighbors of a node i can be also highly connected to the rest of the graph (connections to highly connected nodes contribute more than those of the weakly connected nodes). Hence, the eigen-centrality measures more the global influence of i on the network and verifies:

$$C_i^{eigen}=\sum _{j=1,\dots ,n;j\ne i}{C}_j^{eigen}/\lambda$$

(14)

where $\lambda$ is the greatest eigenvalue of the adjacency matrix A of the graph.

The four centralities above have each their intrinsic interest: (i) the betweenness centrality of a node is related to its global connectivity with all the nodes of the network, (ii) the degree centralities indicate the level of local connectivity (in-, out- or total) of a node with only its nearest nodes, (iii) the closeness centrality measures the relative proximity a node has with the other nodes of the network for a given distance on the interaction graph of the network, and (iv) the eigen-centrality corresponds to the ability for a node to be connected to possibly a few number of nodes, but having secondarily a high connectivity, which indicates for example an important relay for the dissemination of news in social networks, the social networks being also used to model the spread of rumors, political opinions and social diseases, such as obesity and its most important consequence, type II diabetes.

For example, for curing a social disease, an Hospitaller Order of Saint Anthony has been founded in 1095 at La Motte (presently Saint Antoine near Grenoble in France) by Gaston du Dauphiné, whose son suffered from a fungal disease, known in the Middle Ages as Saint Anthony’s fire or ergotism, caused by a transformation of grain like rye into a drug (the ergotamine) in flour provoking convulsions often leading to death. The members of the Saint Anthony Order were specialized in curing patients suffering from ergotism and they spread with about 370 hospitals over the whole Western Europe in the fourteenth century, able to treat about 4000 patients. They observed that the disease propagation was due to social habits like eating same type of bread, then ergotism is considered as the first social disease whose contagion process depends on alimentation behaviour.

Now obesity is considered as the most characteristic social “contagious” disease, and has been defined as a pandemic (i.e., an epidemic having a world wide prevalence) by the WHO . Both nutrition mimicking and social stigmatization explain the dissemination of obesity into a social network [64,65]. Obesity consists in an excessive accumulation of fat in adipose tissue favoring chronic diabetic and cardiovascular diseases. In UK for example, 36% of men and 33% of women are predicted to be obese in 2030 for only 26% of both sexes in 2010 [66,67,68,69] showing the critical character of the present explosion of obesity, especially in young people. For assessing the prevalence of the disease in young population, we have observed social graphs based on the existence or not of friendship between pupils in the 5th and 4h classes of two high schools (corresponding to ages from 11 to 13 years), one in France at Jœuf already published in [5], and a new one in Tunisia at Tunis. The corresponding graphs (Figure 6 Top) and the histograms of the number F of friends of the pupils, show in both cases that the distribution of F is uni-modal for normal weight pupils and bimodal for overweight or obese ones (Figure 6 Middle), because of a double phenomenon in obese pupils: a part of them attracts with a friendship link due to their open personality, and another part repulses. Then, we have done the simulation of a preventive education of central obese nodes following WHO recommendations [70] during the dynamics of obesity spread following the Equation (7), in which we give the tolerance h = 0 (resp. h = 1) to the obese of the first (resp. second) mode. This education has been given in Figure 6 to the obese central nodes for betweenness in the social graph G. The result is the best between the four centralities defined above, but it remains an irreducible percentage of obese in the population of tolerance zero (Figure 6 Bottom right). For that reason, we need to introduce a new notion of centrality based on the notion of entropy.

4. Discussion

4.1. The Notion of Entropy Centrality

Indeed, despite the interesting properties of the various classical centralities (Figure 6, Figure 7 and Figure 8), the therapeutic education targeting the most central nodes leaves in general a not-neglectible percentage of obese after education, then we will discuss now the introduction of a new notion of centrality, called the entropy centrality, taking into account the heterogeneity of the vector having as components the state and the tolerance, into the neighbor set of the node i, and not only the connectivity of the graph:

$$C_i^{entropy}=-\sum _{k=1,\dots ,s}{\nu }_{k}Log{\nu }_k,$$

(15)

where ${\nu }_k$ denotes the kth frequency among s of the histogram of the values of the vector (state, tolerance) observed in the neighborhood $V_i$, set of the nodes linked to the node i.

This new notion of centrality is more useful than the others to detect the good candidates to be educated in order to transform their state to the normal state 0, because they can influence more efficiently a heterogeneous environment to recover the normality. An illustration of this fact is done in the example of the Section 3.3.

By comparing the results obtained thanks to a preventive and therapeutic education leading the overweight or obese nodes to the normal weight, we observe that with only the 21 individuals the most entropic central educated (Figure 8 Bottom right), all the population of pupils is going to the normal weight state, but we need 68 individuals with the total degree, and 85 individuals with the in-degree and eigenvector centralities (Figure 8 Top and Bottom left). Then, the best public health policy against the obesity pandemics consists in using the notion of entropy centrality to select the targets of the therapeutic education [5].

4.2. The Mathematical Problem of Robustness and the Notion of Global Frustration

Each of the examples presented in Section 3 have used the dynamics governed by the deterministic Hopfield transition of the Equation (7) [71]. We have exhibited complex multi-attractor asymptotic behaviors and the central problem to discuss now is the conservation of the attractor organization (namely the number of attractors, the size of their attraction basins, and their nature, i.e., fixed configurations or limit-cycles corresponding to oscillatory states).

This problem of robustness of the attractor landscape is related to the architecture of the interaction graph $G\left(J\right)$ linked to the (discrete or continuous) Jacobian matrix J. The problem of robustness has been often considered since 50 years [41,72] under different names (bifurcation, structural stability, resilience, etc.) and it is still pertinent [73,74].

If we consider the Hopfield rule with a constant absolute value $w>0$ for its non-zero interaction weights, we can study the robustness of the network with respect to the variations of w, by using the two following propositions [51,52]:

Proposition 1.

Let us consider a deterministic Hopfield Boolean network which is a circuit sequentially or synchronously updated with constant absolute value w for its non-zero interaction weights. Then, its dynamics is conservative, keeping constant on the trajectories the Hamiltonian function G defined by:

$$G\left(x\left(t\right)\right)=\sum _{i=1,\dots ,n}{(x_i\left(t\right)-x_i(t-1))}^2/2=\sum _{i=1,\dots ,n}{(g\left(w_{i(i-1)modn}{x}_{i-1}(t-1)\right)-x_i(t-1))}^2/2,$$

(16)

where g is the Heaviside function. $G\left(x\right(t\left)\right)$ is equal to the total discrete kinetic energy, and to the half of the global dynamic frustration $F\left(x\left(t\right)\right)={\sum }_{i=1,\dots ,n}{F}_{i,(i-1)modn}\left(x\left(t\right)\right)$, where $F_{i,(i-1)modn}$ is the local dynamic frustration related to the interaction between the nodes $(i-1)$ and i: $F_{i,(i-1)}=1$, if sign $\left(w_{i(i-1)}\right)$ = 1 and $x_i\left(t\right)\ne x_{i-1}(t-1)$ or sign $\left(w_{i(i-1)}\right)=-1$ and $x_i\left(t\right)=x_{i-1}(t-1)$, else $F_{i,(i-1)}=0$.

The result of the Proposition 1 holds if the network is a circuit for which transition functions are Boolean identity or negation [75] on which we can easily calculate the global frustration F and show that it characterizes the trajectories and remains constant on them (Figure 9). Then, consider now the dynamic entropy E of a Hopfield network [53,76] as a robustness parameter quantifying the capacity the network has to return to its equilibrium measure $\mu$ on the configuration space $\mathsf{\Omega }={\left\{0,1\right\}}^n$ after a perturbation. We can calculate E by using the following equality:

$$E=E_{\mu }-E_{attractor},$$

where $E_{\mu }=-{\sum }_{x\in \mathsf{\Omega }}{\mu }_{x}log{\mu }_x$, ${\mu }_x=\mu \left(\left\{x\right\}\right)$ being the value of the invariant measure $\mu$ at the configuration x, and $E_{attractor}=-{\sum }_{k=1,m\le 2^n}\mu \left(C_k\right)log\mu \left(C_k\right)$, m being the number of the network attractors and $C_k$ the attraction basin of the attractor $A_k$.

$\partial E/\partial w$, which denotes the derivative of the dynamic entropy E with respect to the absolute value of the weight w, can be considered as a new robustness parameter. We have: $\partial E/\partial w=\partial E_{\mu }/\partial w-\partial E_{attracteur}/\partial w$ and we can prove the following result [8,51]:

Proposition 2.

In the parallel updating mode of a Hopfield network, we have:

$$\partial E_{\mu }/\partial w=-wVa{r}_{\mu }\left(F\right),$$

where $Va{r}_{\mu }$ is taken for the invariant Gibbs measure μ defined by: ${\mu }_x=exp({\sum }_{i\in x,j\in y}{w}_{ij}{x}_i{y}_j)/Z$, with $Z={\sum }_{x\in \mathsf{\Omega },y\in \mathsf{\Omega }}exp({\sum }_{i\in x,j\in y}{w}_{ij}{x}_i{y}_j)$.

The Proposition 2 indicates that there is a close relationship between the robustness (or structural stability) with respect to the parameter w and the global dynamic frustration F. This global dynamic frustration F is in general easy to calculate. For the configuration x of the genetic network controlling the flowering of Arabidopsis thaliana [77], there is only two frustrated pairs, hence $F\left(x\right)=2$ (Figure 10).

More generally, all the energy functions introduced in the examples are theoretically calculable and could serve as in physics for quantifying the structural stability. When the number of states in a discrete state space E is too large, or when the integrals on a continuous state space E are difficult to evaluate, then Monte Carlo procedures can provide a good estimation of functions like the dynamic entropy or the variance of the global frustration.

5. Conclusions and Perspectives

We have defined some new tools based on the notion of entropy common to both continuous and discrete dynamical systems, based on thermodynamic concepts. These notions are useful in biological modeling to express in which way precise energy functions related to the various concentrations or population sizes involved in the transition equation of a biological dynamical system can be conserved or dissipated along the trajectories until their ultimate asymptotic behavior, the attractors. Future works could be done in the framework of more general random systems [78,79] in which the notions of entropy, invariant measure and stochastic attractors (called also confiners in [54]) are classical to interpret the data observed in all the fields (neural, genetic and social) considered in the present study and use the interpretations they allow in practical applications (see for example [19,80] for some biomedical applications).