Phase transitions in small-world systems: application to functional brain networks

In the present paper the problem of symmetry breaking in the systems with a small- world property is considered. The obtained results are applied to the description of the functional brain networks. Origin of the entropy of fractal and multifractal small-world systems is discussed. Applying the maximum entropy principle the topology of these networks has been determined. The symmetry of the regular subgroup of a small-world system is described by a discrete subgroup of the Galilean group. The algorithm of determination of this group and transformation properties of the order parameter have been proposed. The integer basis of the irreducible representation is constructed and a free energy functional is introduced. It has been shown that accounting the presence of random connections leads to an integro- differential equation for the order parameter. For q-exponential distributions an equation of motion for the order parameter takes the form of a fractional differential equation. We consider the system that is described by a two-component order parameter and discuss the features of the spatial distribution of solutions.


Introduction
Recently Eguiluz et al. [1] have presented a method of construction of functional brain networks proceeding from the results of functional magnetic resonance imaging measurements in a human. In these experiments, a magnetic resonance activity of certain parts of brain (so-called voxels) is measured at each discrete time step. By we denote a voxel's activity at the instant of time . It was proposed to consider that two voxels are functionally linked if the value of their temporal correlation exceeds a certain positive value independent of the value of their anatomical connection. The correlation coefficient between any pairs of voxels and is calculated as where 2 = 2 ( ) − ( ) 2 , brackets … represents temporal averages and ( ) is the blood oxygenation level dependent signal of the voxel in case of brain scanning data.
The elements of the correlation matrix determine the value of correlations among various parts of the cerebral cortex. Using highly correlated nodes Eguiluz et al. [1] have constructed a network and determined that the degree distribution of the obtained network has the form ~− , where ≈ 2 (figure 1) [1]. It has been also shown that these networks possess a small-world structure, a community structure and are fractals [2].
However it should be stressed that the degree distributions for = 0.6 and = 0.7 are actually described by the -exponentional distribution and for = 0.7, ≈ 1 whereas for = 0.6, > 1 [3]. For = 0.5 the form of the degree distribution of the correlation network is characteristic of stretched exponential probability distributions [3]. Fitting with the help of formula (2) is actually performed by using the maximum likelihood method [5].
In paper [6] the correlation network obtained from functional magnetic resonance imaging measurements is compared with that of the derived from the numerical simulation of 2D Ising model at various temperatures. Near the critical temperature a striking similarity in statistical properties of these two networks are observed and it makes them indistinguishable. The similarity of these networks allows to suppose that a collective dynamics is inherent in a human brain.
On the other hand it also follows that the dynamics of the brain functioning in such systems takes place near the critical point and generates spatio-temporal structures.
To study the dynamics generating spatio-temporal structures in brain we develop the theory of the Landau -Ginzburg type for the systems with a small-world structure.
We first introduce an algorithm to generate small-world networks. Then applying the maximum entropy principle and accounting multifractality of the system we determine a degree distribution for such systems. Using the principle of least action we derive the equation of motion for the order parameter representing spatio-temporal structures near the critical point. We also consider a specific example.

Topology of functional brain network with a small-world property
The algorithm constructing a small-world network was for the first time proposed by Watts and Strogatz [7]. At the initial instant of time there is one-dimensional lattice of nodes with periodic boundary conditions where each link connecting a vertex to one of its nearest neighbours in the clockwise sense is left in place with probability 1 − , and with probability is reconnected to a randomly chosen other vertex. Long range connections are therefore introduced. As a result a network structure with a small-world property emerges. However this network is not scale-free that is it has degree distributions with non-power law forms [8]. A human brain consists of ~10 10 neurons, each of which is connected with other neurons by ~2 • 10 4 links. As is shown in [1] the degree distributions of functional brain networks in the log-log scale have a rectilinear region. Then another construction of a small-world network is a more convenient model for a human brain. We proceed from a closed system of nodes with periodic boundary conditions where each node is connected with the neighbors by links. A new edge is added to this system at each instant of time, one of the ends of this edge being connected with one of the nodes of the regular lattice with probability 1/ while the other end of the edge being connected with that in accordance with the preferential attachment principle / . Thus the rate of the connectivity change of the node (without taking into account the contribution of the initial regular graph edges) is determined by two contributions and represented by the equation Taking into account that at the instant a full network connectivity is = 2 and the change in the total degree of the network at one time step is ∆ = − − 1 = 2 we obtain = 1. Then equation (3) takes the form The solution of this equation has the form where is a constant of integration which can be determined from the condition = 2 . The degree distribution of this network is described by the -exponential distribution in the form [9] Here is a measure of complexity of the system. The value = 1 corresponds to the degree distribution in the form of the Gaussian distribution emerging at large temporal shaping of such system. At ≠ 1 and for large enough the network is characterized by the degree distribution emerging at the earlier steps of forming a network. Distribution (6) describes the results of processing of functional magnetic resonance imaging measurement data = 0.6 [1] pretty well.
The presented algorithm shows that order and disorder are inherent in small-world systems. However it should be taken into consideration that the structure of these networks may be homogeneous, fractal and multifractal.
To derive a distribution function for such systems we use the maximum entropy principle. The famous Boltzmann-Gibbs entropy which is determined as is characteristic for homogeneous structures, where the distribution function is normalized to unity.
We consider the quantity Fractional generalization of this integral gives where is a gamma function. For admissible systems when = in case of a fractal structure we have = ln presenting the Shafee entropy [8]. It is clear that at = 1 we obtain the Boltzmann -Gibbs entropy. A multifractal is a mixture of fractals. Suppose that dimensions of fractals in a multifractal are distributed uniformly that is = 1/( − 1). Then multiplying (9) by and integrating over we obtain Consequently for admissible systems we have representing the Tsallis entropy [7]. If → 1 the Tsallis entropy coincides with that of Boltzmann -Gibbs.
These results can be generalized if we consider a multifractal structure as a mixture of fractal substructures. We introduce the distribution function allowing to introduce an equation for determination of the systems entropy and , is an incomplete gamma function. Considering the standard constraints =1 = 1 and where are observed quantities and applying the maximum entropy principle we determine the distribution function where = 1, … , , and and are the Lagrange multipliers determined from the constraints. Distribution (16) is characteristic of the systems with a skewed degree distribution. Note that distribution function (16) describes the degree distribution of the functional brain network precisely enough in case = 0.5.
Using the constraints = 1 , −1 = , −1 2 = 2 and entropy (11) in the framework of the maximum entropy principle we obtain the degree distribution in the form (6) describing the degree distribution of the functional brain network in case = 0.6, and where = [9]. In case of constraints = 1, = and entropy (11) the maximum entropy principle gives a distribution function in the form where describing the degree distribution of the functional brain network precisely enough in case = 0.7.

Order parameter in functional brain networks
The exact solution of the Ising model has been found only for the one-dimensional and twodimensional cases in the zero external field [12]. So for the description of spatio-temporal structures occurring in the functional brain networks we shall develop a theory of phase transitions of Landau -Ginzburg type [13]. It should be stressed that the theory of phase transitions in graphs is discussed in [14] and analysis of the systems with long-range space interactions and temporal memory is given in [15]. We shall proceed from the construction of a small-world network which is a regular graph with additional shortening links. Under such consideration a regular substructure of a small-world network could possess a symmetry of the discrete subgroup of the Galilean group.
There occurs a problem to determine transformation properties of the order parameter. To solve this problem we first have to determine a number of the components of the order parameter involved into the process under consideration. We proceed from the measured time series and use the Takens method [16,17]. Choosing from a set of experimental data equidistant points, we obtain a set of discrete variables 0 , … , −1 + − 1 , where = 1, … , − 1 and is a temporal shift. It is necessary construct the points of the phase space and calculate a correlation function of the attractor.
where is the Heaviside function and is distance. Further we construct the dependence ln on ln [16,17]. If the value of the slope depending on reaches a plateau above a certain then the system represented by the present temporal sequence must have an attractor. The value which reached saturation should be considered as a dimension of the attractor and the value above which the saturation is observed should be considered as the minimum number of variables necessary for modulation of the behavior corresponding to the current attractor [16,17]. Based on the results of the analysis of the time series we can choose the as a number of the order parameters.
If 0 is an irreducible representation of the discrete subgroup of the Galilean group of a regular subgraph and , = 1, … , , = 1, … , are the order parameters which are transformed in accord with this irreducible representation where is a number of wave vectors in the star * of the irreducible representation 0 and is dimension of a small representation, then the dimension 0 is = × . Thus the choice of a particular irreducible representation could be made from the standard reference books of irreducible representations [18].
It should be stressed that if the wave vector under consideration is expressed irrationally in terms of vectors of the reciprocal lattice then only invariants of the rotation group have corresponding irreducible representations. If the wave vector is expressed rationally in terms of vectors of the reciprocal lattice then the irreducible representation under consideration admits anisotropic invariants.
The phenomenon of hysteresis is inherent in a human brain and thus an irreducible representation may not satisfy the Landau condition that is contain invariants of the third order.
As small-world systems are initially inhomogeneous a free energy functional must have the Lifshitz invariant. So the task is to consider irreducible representations corresponding to the internal points of the Brillouin zone. During the transition into modulated structures the point symmetry elements of the initial phase remain [19]. As the Landau and Lifshitz conditions are violated we deal with phase transitions of the first order.

Derivation of the equation of motion for the order parameter near the critical point in the functional brain network
To derive an equation of motion for the order parameter in the system with a small-world property we determine the free energy functional in the form [13,15]:

(22)
We took into account the presence of random shortening links in the structure. Here is a spatial coordinate, is time and the functions 0 , , ′ , ′ and 1 , , ′ , ′ describe the influence of the small-world property on the critical properties of the system. Integration is performed over the region in two-dimensional space 2 to which , belong. The equation of motion for the order parameter , is derived with the use of the Gateaux derivative of the functional , , determined as: where = is a smooth integrable function. Further it is suitable to introduce the functions: We note that such choice of 0 , , ′ , ′ and 1 , , ′ , ′ allows to separate spatiotemporal derivatives. The dynamic equation for the order parameter is determined from the principle of stationarity , = 0 and for the arbitrary function , has the form: This is an integro-differential equation allowing to obtain an equation of motion for the order parameter for various kernels 0 , , ′ , ′ 0 and 1 , , ′ , ′ .
We suppose that 0 , , ′ , ′ = − ′ 0 , ′ , 1 , , ′ , ′ = − ′ 1 , ′ . In this case time and the spatial coordinate separate. In case of one component order parameter in accord with the symmetry arguments , = 2 2 , + 4 4 , , then from (27) we obtain: Here is a control parameter of the system and is a positive parameter. The form of the function is determined from the condition of invariance. Hence includes the integer basis of invariants of the irreducible representation of the space group of the high symmetry phase. The solutions of (28) are investigated in [20].

Systems described by two-component order parameter
In case of the system described by the two-component order parameter which is transformed according to the irreducible representation corresponding to the internal rational point of the Brillouin zone we have an integer basis consisting of two invariants. In the polar coordinate system these invariants have the form 2 and cos , where and are an amplitude and a phase of the order parameter correspondingly and is a parameter of anisotropy.
In this case a spatial dependence of the order parameter in accord with (27) takes the form: The integral expression in (29) could be considered as averaging over the distribution function where is a renormalized phase of the order parameter. When = 2 the spatial distribution of is determined by the equation (33) In case 0 > 0 2 , the trajectory increases indefinitely.
where am is the Jacobi elliptic function.
The solution corresponding to 0 = 0 2 determines a separatrices. In this case d d = ±2 0 cos 2 , and for the initial condition = 0 = 0 the solution of this equation has the form: For finite motion the solution of fractional equation (30) is determined as is the Mittag-Leffler function and the function + is the Euler gamma function. When = 2 we have 2.2 − 2 = sin and the solution of the linear fractional equation has the form = sin . When 1 < < 2, we have a continual number of differential equations. The oscillating function with a decreasing amplitude is a solution of (30).

Conclusion
Order and disorder are inherent in small-world systems and consequently in functional brain networks. We propose the algorithm constructing a network equivalent to functional brain networks. We suppose that a regular substructure of a small-world network could be described by a discrete subgroup of the Galilean group. In such approach random connections generate the medium in which the order parameter moves.
We have presented the equations allowing to determine an entropy for fractal and multifractal small-world networks and applying the maximum entropy principle we have derived characteristic degree distribution functions. The obtained various distribution functions in the log-log scale have a rectilinear region. Using this property with the help of the principle of the least action we have derived an equation of motion for the order parameter in the form of a fractional differential equation.
We have presented the algorithm for determination of the transformation properties of the order parameter. The Rapp diagram is constructed from the analysis of the measured time series. The number of the components of the order parameter is determined from this diagram. On the other hand the number of the components of the order parameter coincides with the dimension of the irreducible representation equal to that of the irreducible representation of a small group multiplied by a number of wave vectors in the star corresponding to the wave vector under consideration. These data can be found in literature [13].
We show that in case of the functional brain networks only internal points of the Brillouin zone should be considered. Therefore we deal with the phase transitions without the loss of the point symmetry elements [14]. We show that in case of two-component order parameter the spatio distribution of the order parameter is determined by the one-dimensional sine-Gordon fractional differential equation. The change of the fractional dimension due to the fractal dimension of the structure changes the medium in which the order parameter moves. It could be shown that if the fractional dimension of the order parameter equation of motion changes from 2 to 1 the solution different from zero becomes identically zero. Thus our results discover the principles of a human brain functioning near the threshold of criticality.