Phase synchronization and criticality in a network of neural mass models

Synchronization has an important role in neural networks, and their dynamics are mostly accompanied by cognitive activities such as memory, learning, and perception. These activities arise from the collective neural behaviors and are not totally understood yet. This paper aims to investigate a cortical model from this perspective. In this paper, we investigated a network of neural populations in a way the dynamics of each node corresponded to the Jansen-Rit neural mass model. First, we put this dynamic on a single mass of four different input levels. Then, we considered a Watts-Strogatz network of Jansen-Rit oscillators. We observed an epileptic activity in the weak input level. The network to change various parameters is considered. The detailed results including the mean times series, phase spaces and power spectrum revealed a wide range of different behaviors such as epilepsy, unimpaired, and a transition between synchrony and asynchrony states. Since the critical state is a dynamic candidate for healthy brains, we considered some measures of criticality and investigated them on phase transition points. We showed that the criticality hypothesis is not all or nothing theory. It means that due to the nature of speciﬁc indicators selected for studying the criticality, the phase transition point can be a critical point or not. Indeed, some markers of criticality can exist in phase transition points, and others may not. As a result, we do not claim that neural models show criticality or not, but we can only admit that they have a percentage of criticality.


Introduction
Synchronization is a process that the elements of the system, coordinate one dynamical property by coupling or external force. The simplest of synchronization is full synchronization proposed by Fujisaka and Yamada in 1983 1 . In 1996, Rosenblum was the first one who introduced the concept of phase synchronization 2 . Due to this definition, the phase difference will be fixed (smaller than 2π), but amplitudes not necessarily the same in this definition. At least two elements exist in the synchronization process. To investigate network synchrony, we first assume that N oscillators interact with each other. Each oscillator state shows a fixed-point attractor (steady-state), a limit cycle (a closed cycle that shows periodic oscillations), or chaos (a complex orbit that determines aperiodic oscillations). Measurement of synchronization in neural signals is a procedure to estimate oscillatory activities [3][4][5] . The rhythmic neural interaction can measure by different linear (e.g., cross-correlation) or nonlinear (e.g., mutual information) methods that analyze signals in time and frequency spaces 6 . Synchronization measures have been wildly interested in investigating neural dynamics [7][8][9][10] . Especially, synchronization phenomena have a substantial role in determining normal 11,12 and abnormal 13 brain function and also in a deeper understanding of information processing [14][15][16] . Moreover, synchronization influences on prediction and detection of some disorders such as epilepsy 17 , Alzheimer 18 , autism 19 , and schizophrenia 20 . About 1% of people suffer from epilepsy 21 . Many patients can control their disease by taking medicine, but approximately 30% of them have drug-resistant epilepsy. The dynamics of epilepsy disorder is very complex, and it arises from high synchrony of neuron activities 22 . The significant changes in the brain dynamics during epilepsy can be fatal, and in certain conditions, lead to loss of consciousness, body tremors, and even death. To enhance knowledge of the mechanism of epilepsy is imperative. The biophysical modeling of the brain is an active area for neuroscientists. In this approach, they simulate a complex system by mathematical and physical tools to predict the influence of different biological factors [23][24][25] . This attitude causes a profitable framework in neural dynamics modeling with large-scale approaches [26][27][28][29] . Some cognitive activities such as memory, learning, and perception arise from a mass activity of neurons. There are different scales to study the brain function between a small patch of cortex and a large-scale system. Coupling a collection of neural mass models into mesoscale circuits can bridge between these scales. The population of neurons coupled together to form a neural mass and they can model different phenomena and generate a variety of their dynamics behaviors.

Graph Theory
Modeling the brain as a complex network is a powerful mathematical tool to understand the structural and functional of the brain architecture. Structural brain networks can be shown as graphs that their nodes (vertices) are related to neural elements (neurons or brain regions) and linked by physical connections (edges). Using adjacency matrix (M) is a simple method to display a matrix where M i j = 1, if there is a connection between node i and node j and Mi j = 0 otherwise. In this work, neural populations are considered as nodes, and edges are interpopulation connections. The clustering coefficient computes how connected a node's neighbors join together to make a cluster. This feature in a network is defined as the average clustering coefficient of nodes. The shortest path length is the least number of edges between nodes and the average of this property shows the average node to node distance in a graph. This concept represents how rapidly information can be transformed through the network. Regular, random, and small-world networks are three important network models which have largely been studied so far. In regular networks, each node has exactly the same number of connections and the clustering coefficient is high in these networks. The connections between each node in random topology follow a normal degree distribution with a low level of clustering. The average path length is short (long) in regular (random) systems. Regular, random, and small-world networks are three important network models which have largely been studied so far. In regular networks, each node has exactly the same number of connections and the clustering coefficient is high in these networks. The connections between each node in random topology follow a normal degree distribution with a low level of clustering. The average path length is short (long) in regular (random) systems. Evidence suggested that the real-world networks have small-world properties including two independent structural features, namely, short average path length (long-range connections) and high clustering coefficient. The Watts-Strogatz (WS) models were generated as the simplest networks that have the small-world properties.

Jansen-Rit model
In this paper, neural masses are defined as columns. The model consists of a main population of pyramidal neurons that receive inputs from three resources: excitatory and inhibitory interneurons feedback in the same column and external input from other columns. Figure 1 is a schematic representation of this model. We consider a network with an interactive population, representing a patch of the cortex. Each node indicates a neural mass with Jansen-Rit dynamics. Nodes can interact with each other and generate different behaviors. First, we represent this dynamic per node k, k =1,. . . , N. The equations are written by three second-order differential equations and can rewrite with six first-order differential equations as follow 54 : Where (y0, y3), (y1, y4) and (y2, y5) are activities of pyramidal, excitatory and inhibitory ensembles respectively. S is a sigmoid function, transforms the average membrane potential of neurons to mean firing rate of action potentials, and identified as follows: Average numbers of synapses between excitatory neural populations 1 * J, 0.8 * J J3, J4 Average numbers of synapses between inhibitory neural populations 0.25 * J J Average numbers of synapses between neural populations 135 v max Maximum firing rate 5 Hz v0 Potential at half of maximum firing rate 6 mV r Slope of sigmoid function at v0 0.56mV −1 (y1 − y2). Also J i = Jα i , for i =1, . . . , 4. Now, the Jansen-Rit's dynamical equations in a network with N nodes, i=1,. . . , N, are as follows 41 : Where M i j shows the adjacency matrix of network and α represents coupling coefficients between nodes.

Pearson cross-correlation
Neural interactions can be measured by multiple criteria, each with its advantages and disadvantages 45 . The simplest measure of non-directed interactions between random variables is Pearson cross-correlation (Pcc) which measures the linear relationship between each pair of random variables. This coefficient removes the temporal structure, so time series are considered as generalizations of random variables. Normalized Pcc coefficient in one-dimensional between two time-series x and y is defined as follow: Where N shows the length of the signal and < x > is the mean of time series x. The range value of this coefficient varies between -1 and 1. The quantity of 1 (-1) represents a perfect linear positive (negative) correlation and 0 means no correlation between two series (uncorrelated state). Note that the matrix constructed by this definition is symmetric and elements are one in the main diagonal.

Detrended fluctuation analysis (DFA)
DFA method specifies the self-affinity of a signal and demonstrates long-range correlation in time series, which was first introduced by Peng et al in 1994 46 . Illustration of this procedure consists of the following five steps: • Step 1: We consider time-series x with a length of N and calculate the cumulative sum of where t ∈ N.
• Step 2: Y t is divided into N s = (N/s) segments with the same size s. It is possible that (N/s) not to be an integer, so we do this division from the end of signals. Therefore, we have 2N s segments. Then, the least-squares fit is calculated by minimizing the squared errors in each segment.

4/20
• Step 3: In each segment, the root mean square (RMS) for Y is calculated as follow: which Y local shows the local trend.
• Step 4: We repeat step 3 for different interval sizes ranges between 4 and N/10 to find the relation between the scaling exponent and the data fluctuation. • Step 5: We plot F(s) consistent with s in a logarithm scale. In power-law relation, F(s) ∼ s a which a, fluctuation or scaling exponent, is determined as the slope of a straight line fit. 0.5 < a < 1 represents long-range correlation.
We used Nolds package in Python with no overlapping windows to perform this technique (https://pypi.org/project/nolds/).

Simulation
In this simulation, we considered a small-world network (Watts-Strogatz network with rewire probability equal to 0.2) with 80 oscillators as nodes. Each node is connected to 20 neighbors, 10 on each side. We applied Jansen-Rit dynamics in all nodes and used the stochastic Runge-Kutta method to simulate the system 67 . The time step was set to 10 −4 s. The coupling coefficient between neurons in each mass was considered 135 (to show alpha rhythm in each single node). Each node received external inputs from other brain regions. First, we supposed that the coupling coefficient between nodes was equal to zero and then increased it. External input interval was [120, 320] which we divided into four segments. Based on the permitted input interval, we analyzed the dynamics of the network for four different inputs named small, intermediate, strong, and ultra-strong input. For each of them the input was Gaussian random noise with mean of 145 , 195, 245, and 295 respectively and sigma = 3.25 for all of them. The length of this simulation was 200 seconds and ran ten times on each coupling coefficient. Elimination of the first 5 seconds of simulation, gave us the certainty that the system achieved to its equilibrium.

Results and discussion
Synchronization phase transition First, we applied the Jansen-Rit model on a single node. The power spectrum of the output signal for four levels of input is shown in figure 2. Increasing input increases the dominant frequency, but all of them are still in the alpha band. Then, we  consider a network with an interactive population, representing a patch of the cortex. Each node shows a neural mass with Jansen-Rit dynamics. There are a couple of factors that can influence the network's dynamics. One of the most important of this is coupling strength. Nodes can interact with each other and generate different dynamics behaviors. Changing the coupling coefficient between nodes and analysis of network dynamics has been interested recently 35,38,39 .

5/20
At the start, we assumed that each node received the external input from a weak range. Cross-correlation is a popular method to investigate and compare time series. By considering this time series of each node in the model and calculating Pcc between them, a square matrix was created. The ij'th element of the matrix shows the Pcc of node i and j dynamics. The mean of the Pcc matrix for weak input is shown in figure 3. Vertical red bars represent dispersions for ten runs (standard deviations). No coupling leads to zero cross-correlation. It can be seen in figure 3 increasing coupling up to 1.5 causes more correlation. In addition, the mean signal amplitude gets larger (first column in figure 4). During this range of coupling strength, the dominant frequency is still in the alpha band (α = 0 in figure 5(A)). At a coupling coefficient equal to 1.5, the network's behavior abruptly has been changed, and as it can be seen in figure 6(C), the large amplitude oscillations appear. In each run, switching from high to low amplitude oscillations occurs at different times. In addition, its phase space is exhibited in figure 6(D). It is visible that in this coupling coefficient, the phase space is composed of 2 limit cycles. The first of which in respect of high amplitude oscillations show rhythmic spiking activity, and the second corresponding to a shorter amplitude signal. This transition between high and low amplitude signals is one of the most indicators of brain disorders such as epilepsy [68][69][70] . High amplitude signals are rhythmic spiking activities represented in theta band frequency and show epileptic activity that the frequencies of recorded epileptic seizure confirm this 68 . As well as, the dominant frequency of the short amplitude regime is alpha. Indeed, this behavior is mixed epileptic-alpha activity. It means that the ictal activity occurs 58 and it gradually slows down in a non-specific time during the whole simulation due to the stochastic noise and different initial conditions in each run. As a matter of fact, a spontaneous transition happens between two types of behaviors. Results indicated that in the [1.5, 3] regime of coupling strength, the output signal for all repetitions shows only two types of behavior: high amplitude oscillations that either change to low amplitude oscillations or not ( figure 6). Notably, the length of this simulation was 200 seconds. Indeed, epileptic activity is seen in this area that can be continued up to 200 seconds or change to a low amplitude oscillatory state. Consequently, we call this area an epileptic regime, and the coupling coefficient equal to 1.5 (3) is the starting (ending) point of epilepsy.
Interestingly, that bistability and transition between periodic areas with high and low amplitudes are the two most features in the investigation of epilepsy 69 .
Fast Fourier Transforms (FFTs) for the values of 1.5, 2, 2.5, and 3 show a significant change in their appearance ( figure  5(A)), which confirms a disorder in network function 69,71,72 . Figure 5(B) shows the dominant frequency according to different coupling coefficients. The color of each cell represents the dominant frequency of an independent in-silico experiment (run). The vertical axis is corresponding to the coupling strength value of each run in its raw. There is a jump phenomenon in frequency in the [1.5-3] regime that is a key feature to realize the epileptic seizure dynamics 69,73,74 . Exclusive of this phenomenon, it is clear that the coupling coefficient and dominant frequency have an inversion relation. Also, presencing bistability in frequency at α=3 is an interesting result. Based on figure 3, a small increase in coupling coefficient more than 1.5 causes a dramatic rise in the mean of the Pcc, and the network has been brought into the maximum value of synchronization. Besides, it can be viewed that there is a big dispersion of the mean of the Pcc matrix in bifurcation points (1.5 and 3). Actually, in the bifurcation points, the network's dynamic changes dramatically. In coupling coefficient equal to 3, the mean of the Pcc has a great variety, and the variance of 6/20  D, F). The x-axis and y-axis show the output signal and its derivation, respectively. By increasing the coupling coefficient, the mean signal amplitude gets larger. The inset of (A), (C) and (E) shows a specific segment of signal for 3 seconds. them is high. There are more fluctuations of the mean of the Pcc in 3 rather than other coupling coefficients. Cross-correlation matrices in the value of 3 are shown in figure 7. There are different patterns in these matrices. The red (blue) color shows full synchronization (anti-synchronization). Also, in this figure, local and global synchronization can be seen. This issue is determined with red and blue masses, respectively. Local synchronization means that some ensembles of neurons have the same behavior, and some of them act completely differently, i.e., synchronized clusters. Simply put, full red Pearson cross-correlation matrices mean that the network shows a global synchrony state, and the simultaneous presence of blue and red masses in a matrix is a symbol of in-phase and anti-phase local synchronization. This different pattern could be reminiscent of different synchronization patterns of brain networks during different tasks 38 . Again, we followed the rise in coupling strength. It is visible that the decrease in the mean of the Pcc is continued up to 3.5. Figure 8 displays that the mean signal output goes back to the normal state and put in the alpha band frequency (α = 3.5 in figure 5 (A)). More increase in α from 3.5 to 5 causes a growth in the mean of the Pcc up to reach maximum value. If we increase α further, this synchrony behavior switches to a fixed-point state, i.e., each node leaves the limit cycle and collapses to a resting state (fixed-point). Mean signals and their phase space for one repeat are shown in figures 8 and 9, respectively. Moreover, their power spectrum for α = 3.5, 4, 4.5, and 5 in figure 5(A) shows a switch from disorder to the normal state, and the dominant frequency of the network goes into an alpha rhythm. So far, the amount of external input of each node was assumed from a weak level. Here, we consider the input value in the other three ranges (medium, strong, and ultra-strong) and report gathered results. The plot of the mean of the Pcc matrix for four inputs versus coupling strength has shown in figure 10(A). Moreover, for better understanding, figure 10(B) is a heat-map of the mean of the Pcc matrix for ten independent runs in each of them, where the x-axis is coupling coefficients, and the y-axis is input signals. For all different intensities of input, the correlation between units for the zero couplings is too low. Because in this case, oscillators can not have a significant effect on each other. We expected that the weak input level went forward as specified expected behavior trend in [1][2][3] in figure 10, however the dynamics of the network due to set parameters seem very complex and unpredictable. It is obvious that lower intensity input needs a higher coupling coefficient to reach its Pcc maximum. Each unit receives some inputs from its neighbors plus the stochastic amount chosen from an interval. Strengthening coupling between nodes leads to receive more from neighbors, and as a result, less stochastic value is needed. Indeed, the input is a fundamental element, and its variations can impress the amount of correlation and the rate of the upward trend to reach the highest degree of this. Interestingly, though plots in each intensity of input have been shifted to the left side, the dispersions of the mean of the Pcc matrix at the α=1.5, 3 are still high. According to figure 11, one impressive event is viewed: a jump frequency in the weak level input that is related to the epileptic regime and has been confirmed in 69 . Furthermore, in each input level (except in weak), increasing coupling coefficients result in the reduction of the dominant frequency but remain in the alpha band frequency.

Power law and criticality
It claimed that the normal brain acts in a critical regime [47][48][49][50] . It is an interesting and challenging question that how the critical dynamics can be detected in neural models. There are some measures of criticality that explain as follows. The expression of the criticality hypothesis is that the brain may be taken on the edge of phase transition. We considered synchronization phase transition in our research, and so we study α=1.5, 3 candidate points that can be checked criticality on them. It is shown in 75 that a different range of patterns can be generated by being on the edge of synchronization. Last two seconds of the spatiotemporal matrices in α = 1, 1.5, 2, 2.5, 3, and 3.5 are seen in figure 12. In α = 1, the fluctuations are random, and the value of synchronization is low. Conversely, in alpha = 2, the activity of the network is extremely ordered. This pattern is related to full synchronization that does not show a healthy state. Between these points, in α = 1.5, the network activity is among high and low synchrony. Similarly, in α = 2.5 and 3.5, activities are highly ordered and disordered, respectively, and in α = 1.5 and 3, the network shows a pattern between ordered and disordered states. So, it is possible to say these points (1.5 and 3) are critical. Another marker of criticality is the presence of high fluctuations. In critical points, the perturbations Figure 12. The spatiotemporal matrices of the mean signal for different values of coupling coefficient. In α = 1.5 (3), there is a phase transition from disordered (ordered) to ordered (disordered) states, and the system in these coupling strengths is being on the edge of phase synchronization. grow in magnitude, and different behavior occurs. Because the system stability is weak, it is expected to see events in each scale (micro or macro) that it is one of the most remarkable properties in criticality 76,77 . Broadly speaking, the border between scales is not clear. Based on figure 3, α = 1.5 and 3 has more dispersion than other coupling strength. An extensive range of synchronization has just arisen in α = 3. Different patterns in the correlation matrices are seen in figure 7. In alpha=1.5, the dispersion of synchronization is approximately high, but the variety of schema fails here. Hence, alpha = 1.5 may not be considered a critical point and does not have a plausible variation of observations. Scale-invariant of neural avalanches is the most hallmark of criticality 78 . In phase transition points, the distinction between microscopic and macroscopic dynamics will be faded, i.e., neural avalanches start on the microscopic level and become so large until they cover overall dynamics and encompass the macroscopic scale. This event has corresponded to the powerlaw distribution of size and duration of avalanches 79 . We analyzed the squared Hilbert amplitude, which is a measure of instantaneous power. We set a threshold (mean of activity (θ )), and we defined a neural avalanche size as the fluctuations of the integrated network mean activity when exceeds a threshold. Also, the duration of neural avalanches, i.e., the time between the start and end of a neural avalanche, is the length of time of the fluctuations stay above θ . The distribution of size avalanches in the bifurcation points for one arbitrary realization in the weak input level is illustrated in figure 13. Although both coupling strengths can generate unusual events during the avalanche, the power-law behavior has not been seen in them. So based on this characteristic, these points are not critical. The presence of LRTCs in the amplitude of neural oscillations supports the critical Figure 13. The distribution function of the size of avalanche in the weak input level at α = 1.5, 3. Neither of them shows power-law behavior.
hypothesis. Indeed, long-range temporal correlations are a pivotal feature of criticality 51,52 . The temporal correlation structures of the signal are investigated by the mean auto-correlation or DFA method. A signal can show LRTCs if its auto-correlation decays as a power law (with an exponent between -1 and 0). Generally, auto-correlation functions are very noisy in their tail, and so, the exponent estimation is very complicated. DFA is a proper technique that overcomes these problems 52 . This technique in neural mass models is used recently in 35 . According to the previous section, it seems that one phase transition occurs in the weak input level around α = 1.5, 3. We checked the presence of LRTCs at these points in an arbitrary repetition. Figure 14 shows the DFA applied in absolute signal Hilbert transform with no overlapping. The linear fit is not suitable for data (on the report of 53 spline fit is the best fit model). So, LRTCs do not exist at these points. Similarly, LRTCs have not been exhibited in other coupling coefficients, and consequently, the network with defined parameters does not show power-law behavior and hence criticality due to this feature.

Conclusion
Our results showed that in a weak input level, the dominant frequency of a single unit and a network of units might be different. In 55 Spiegler. et al. demonstrated that the dominant rhythm could be scaled by the ratio between the inhibitory and excitatory time constants. Although the alpha rhythm was prominent in a single node, the network made by single nodes not necessarily showed this rhythm. This consequence is an emergent property in a complex system 80 . i.e., each isolated node represents an alpha rhythm that can be changed when these units interact with each other in a connected network. We explained elaborately that the input level and coupling coefficient were crucial parameters. Surprisingly, in the weak range of input, some interesting events occurred, e.g., an epileptic activity appeared that related to the first peak in the mean of the Pcc matrix and was unpredictable due to the plot's behavior. α = 1.5, 3 were set as bifurcation points, and [1.5-3] was called "epileptic area". In α = 1.5, GSWDs (generalized spike-wave discharges) with approximately 4 Hz frequency band, were appeared. Commonly, GSWDs have arisen after paroxysmal oscillations in the corticothalamic networks but actually, their process is not clear yet 81,82 . Also, we observed that in an epileptic brain, increasing the coupling coefficient stimulates the network and nonlinear dynamics causes a jump between small and large amplitude oscillations. An important point to note is that different frequency bands have a remarkable role in epilepsy 83 . In α = 1.5, we discovered a transition from preictal to ictal state that is recognized by rhythmic spikes in a theta-band frequency and triggers the seizure activity. It is noteworthy the frequency jump is a vital issue to comprehension epilepsy 69,84 . Until now, all of these results are received from a weak input level. Although a high variance of the mean of the Pcc had shown in α = 1.5 and 3 in intermediate, strong, and ultra-strong input levels, coupling strength and input levels can not produce all behaviors in the weak input level, and so no bifurcation and phase transition exists in these cases. In the final section, we presented some markers of criticality and investigated them on bifurcation points. Results showed a number of the criticality indicators are held in these points and others not. i.e., at the α=1.5 and 3, based on some measures system shows signatures of criticality, and based on some other measures, it does not. Indeed, there is a difference in critical conception in the brain and the statistical mechanics. A critical point in the brain can extend to a critical area, Griffith phase 85 , and criticality is not all or nothing attitude. Moreover, our research has not proven the accept or reject the critical hypothesis. It seems that the Jansen-Rits model has a percentage of criticality theory. As a summary, the randomness of input, initial conditions, and coupling strengths in a network can generate different complex behaviors. We also should mention these behaviors are not just grounds to triggers the seizure activity because understanding the epilepsy dynamics can not be achieved easily. It seems that the delay between the neural populations in networks can produce complex dynamics, and in this study is assumed to be zero. The findings have been reported by 35,86 were considered the delay in other neural mass models. Notice, this model is analyzed in the alpha frequency band. It means that the coupling between excitatory and inhibitory mass in each node was set to 135. Future reports can focus on the epilepsy state, and their results can be effective in the detection and treatment of this disease. Changing the standard deviation of stochastic input may have a significant effect that was discarded here, and maybe some interesting results will be achieved by considering this. Criticality is an asserted assumption in a normal brain. The results in this paper show that the criticality hypothesis is not all or nothing theory. The Jansen-Rit dynamic as a neural mass model with assumed parameters does not have all critical features, and we can only admit that they have a percentage of criticality. We are left with the question, whether the Jansen-Rit model can be modified to show all markers of critical behavior and, if it is possible, adjusting which parameter (noise, delay, coupling between inhibitory ensembles, etc.) do it.