Emergence of Functional Cortical Patterns of neurons characterize the self-organizing way to cognition in brain

Synaptic plasticity and neuron cross-talk are some of the important key mechanisms underlying formation of dynamic clusters of active neurons. In our proposed neuron activity pattern model assorted range of coupled interactions generated the ordered pattern dynamics, mimicking the functional clusters of neuronal activity. The temporal correlates of neurons at near critical synaptic strength features emergence of functional cortical patterns (FCPs) symbolizing task-specific neuronal activation. We investigated the near critical dynamics of long, critical and short-ranged interactions and lead to interpretation that criticality is a range and not just a point. Our model has deciphered the mechanism of emergence of a complex functional brain network. The range of interaction among neurons play an important role in defining the functional severity of the neuronal circuitry. We have found that the coupling range as a function of synaptic strength could make connections in an augmenting fashion from short to long range, as the synaptic strength increases in the population of neurons. Thus, could give insight towards the intensity of cognitive behaviour, as a symbol of multiple stimulus attempts.

The phenomenon of global features of complex cognitive functions in brain arise from the organization of basic functions systematically coordinated in localized brain regions [1]. These two aspects of brain organization are reflected in brain cortical functions and evolve continuously through internal and external perturbations. Hence, human cognition can be thought of as a water surface where patterns in the form of ripples dwell on its surface upon receiving perturbations from outside. The fluidity of consciousness is reflected through the highly dynamic behaviour of brain activity [2]. We can merely access those spatio-temporal patterns of brain activity through various electro-physiological and neuro-imaging techniques in order to understand the underlying dynamics. Electroencephalography (EEG) is one such technique to measure electrical brain activity with high temporal resolution. The functional spatio-temporal clusters of neurons in the cortical layers of the brain is a consequence of various biological factors affecting concentration of neurotransmitter, transfusion of ions through ion channels and voltage generated across the membrane of pre and post-synaptic neurons [3,4]. These biological factors determine the strength of the diffused signal among neurons and neuronal organization. Such recurring signals strengthens the circuitry of the neuron cluster allowing it to learn and memorize a specific function in the brain and rewiring among these complex functional patterns lead to varied cognitive abilities of the brain [5][6][7][8][9]. But, how the synchronous neuronal activation and generation of functional patterns leads to execution of a cognitive task is a riddle to solve. The basis of understanding the complex cognitive functions would begin through studying the network dynamics of the interacting local and global neurons. In shaping the functional network topology of the brain, the interaction range of the cortical functional regions in the brain and neurons within them plays a significant role. And it would be interesting to know that how this range is determined for information processing and shaping neural circuitry among the neurons. In our foregoing work, we proposed a neuron activity pattern model, in which the interaction range was expressed as a function of coupling strength (J) in the modified long-range spin Ising model and categorized interactions into long, critical and short range [10]. It has been showed in [10] that each interaction/coupling range has their unique coupling strength for generating spin clusters at the critical point of phase transition. However, their overall dynamical behaviour has followed universal scaling laws. Our anticipation led to the conjecture that cortical neurons might follow a second-order phase transition near criticality to generate patterns of active neurons in response to a signal. These spin clusters may correspond to mimic neuron activity patterns generated while performing a specific task in the brain. However, the emergence of these functional clusters and their organization are still an open question.
We address important features of spin cluster formation, their organization and properties at below, near and above their critical parameter (Temperature) strengths for each coupling range to generate dynamic patterns of ordered activity which could mimic with the organization of cortical functions of brain. To compare with our model data, we used the EEG data obtained from the UCI database which are recorded while performing a specific visual task by human subjects [11]. We study the topological properties of the network dynamics of such patterns and correlated with the functional networks of the human brain cortex data. In this comprehensive work, we have validated our insights of brain functionality by characterizing the topological properties of self-organized functional cortical patterns embedded in the time-series data of task-specific activity.

Results and Discussion
The Ising model is a simple physical model that could explain the ferromagnetic behaviour with domain coarsening of the two-spin lattice system through inter-molecular interactions imprinted by some range and concurrent computation of the Hamiltonian energy as a driving factor for system evolution. The neuron activity patterns could be imprinted in a connectivity matrix based on correlation in the activity of neurons. We have applied different approaches to build the connection matrix out of the model and empirical time-series data as mentioned in the work flow-chart as shown in Fig.1A. We tag the functional cluster of synchronous activity in the model generated network as functional cortical pattern, FCP, and corresponding to this visual task based EEG time-series has been transformed into functional brain network, FN, using the visibility graph approach [12]. We acquire significant congruence of the characteristics of our model and EEG data in their topological properties and functional characterization through the network theory and multi-fractal approach [9,13]. The model simulations at various interaction ranges and temperature strengths substantiates the emergence of functional patterns of activity at near-critical strength only as depicted in the recurrence plots in Fig.1B, will be explained later in detail. We first need to delineate the details of the model simulation in the next sub-section, succeeded by the deliberation over results.
Correlated neuron activity could optimize brain function Thinking of neural dynamics from a statistical mechanics window is bit challenging but is ruling out since Ising model was used to explain the dynamical avalanches of neuron activity in the brain [2,6]. In comparison with the Ising model a local neural network is not an equilibrium system as neurons constantly receives inputs from other neurons but the near-critical point of phase transition i.e. far from equilibrium could explain the dynamic equilibrium in complex systems. Also, the self-organizing dynamical systems expounds that brain works at a transition point when it receives the optimum signal strength [14]. Still, there is thirst for digging more into the complex dynamics of brain functionality. How brain does the transition into a dynamic functional network that spans the cortical layers of brain acquiring non-random synaptic connectivity for performing a cognitive task [15].
The Ising model is one of the primitive models known to study phase ordering dynamics in a random distribution of spins under the influence of a critical parameter i.e. temperature (T). It had been used profoundly to explain the functional neuronal activity across functional brain networks [16,17]. This makes us strongly presume that phase transition state is the critical state of brain functional dynamics at which cognitive actions take place. In our proposed long-ranged spin ising model (LSM), varied coupling strengths as a function of the coupling range has defined long, critical and short ranged interactions among spins using the Monte carlo Glauber kinetics with non-conserved magnetization [10]. We have considered random allocation of two-spin configuration in a 2D-square lattice with periodic boundary condition in order to avoid the edge effect. We can think of cortical neuronal circuits coarse-grained into the firing and non-firing states of neurons (or group of neurons) represented by the two states of a spin s, s = +1 for firing and s = −1 for rest or non-firing neurons. The total energy of the system representing various internal fluctuations and diffusive processes in neurons that results in their activity and influences the activity of surrounding neurons is given by the following Hamiltonian to model long as well as short range interactions, Where s i and s j denotes the activity of neuron at site i and a neighboring site j and J represents the coupling strength among neurons defined as a function of distance between two ith and jth neurons given by, r ij = |r i − r j | which is confined by n, the coupling range in the system. Then one can model J by power law bahaviour given by, J(r ij , n) = J r n ij [18], where, d is the dimension of the system, and n can be expressed as n = d + σ = 2 + σ [19,20] which can qualify for short range for σ > 2, whereas for long range for 0 < σ < 2 [21,22]. Further, J > 0 always (ferromagnetic case) and n > 0 [18]. We have taken a 2D square lattice of spins with side L and N = L × L being the total number of spin sites. To define coupling range we have assumed centre of the 2D square as the reference For Long (n = 3), Critical (n = 4) and Short (n = 6) coupling ranges at synaptic strengths below-critical (T = 1), near-critical (T ∼ Tc) and above-critical (T > Tc) strengths ( as shown in black, red and green respectively), in a network of 4096 nodes. Note :-There is formation of communities of active neurons at near criticality for all the coupling range interactions.
origin to delimit the maximum distance for the spins to interact. The incircle of radius r = L/2 is the approximated area of the spin distribution which we have considered to calculate the potential energy, U, of the system. In order to capture thermodynamical parameters of the system, one can define the following potential function U (n) as a slowly decay potential function [2,18], where n = 2 + σ for 2D system, scaling J → J/N in the limit N → ∞ and using Euler-McLaurin sum formula [23], we have [10], Critical parameters of the 2D system can be calculated from the force derived from the potential in equation (2), at the critical point where singularity arises in the solution of the system [24]. Following this technique one can estimate critical temperature from the numerically derived closed form approximated equation [24] given below, where, f 2 = <i,j> J(r ij , n) 2 , which can be calculated for the 2D system by following the procedure mentioned above i.e. scaling J → J/N in the limit N → ∞. Now using Euler-McLaurin sum formula [23], f 2 is given by, Now using equations (3), (2) and (4) we could able to arrive at expression for finite critical temperature given below, However, there are debates on the estimation of T c dependence on J in 2D system. Even though the relation between T c and J was obtained analytically in 2D system as sinh(2J/k B T c ) ≈ 1 [25], numerically calculated J/k B T c using Monte carlo techniques in such systems within finite size scaling formalism was found to be varying in the range J/k B T c → [2.269 − 2.29] [26]. The reason for variation in the critical parameter value could be due to finite size (N ) of the system [26], and T c depends on size N except for short range interaction (n > 4). The critical temperature T c can be defined in terms of J/k B , however, through our model simulations of random spin lattices we have got the float values of T c for various long and short interaction ranges (see Fig.1 of ref. [10]), signified by the sudden drop in the total magnetization (order parameter) of the system. This shows that J is the coupling strength which is a factor that controls the temperature T and thus also the critical temperature T c , which could signify the overall synaptic strength that determines the intensity of coupling among the population of neurons. We have applied global synaptic strength to the system of N neurons. This could mimic the synchronous simultaneous activation of neurons sharing common synaptic strength. The critical temperature T c , which signifies the state near the phase transition, could able to trigger optimum strength J to coordinate short and long range neurons in the brain for the generation of functional clusters of neurons at local and global levels. This T c could also mimic the optimal perturbation (stimulus) imparted to the system i.e. relayed and associated in the form of synaptic strength in the brain. From the above equations (5), we could able to understand that for fixed systems size N both short range and long range interaction of neurons contributed to J lowers in T c allowing to work the neurons at modified phase. On the other hand, for fixed n, T c increases as N increases both in critical and long range regime We have taken a 2D 64 × 64 lattice system of 4096 spin sites depicting either firing or non-firing neuron state, allocated randomly. Our approach examines the inter-neuronal connections for long, short and critical coupling ranges at different global synaptic strengths (T ) as categorized into three stages, below criticality, near criticality and above criticality. The below criticality is defined as the temperature below the phase transition temperature but above the mean field temperature i.e. 0 < T < T c . We have considered T = 1 J kB as below critical temperature for all coupling ranges. However, the near critical transition temperature, T c , is different for the long (T c ∼ 2.9 J kB ), critical (T c ∼ 1.9 J kB ) and short (T c ∼ 1.4 J kB ) ranged interactions, based on their respective dynamics (see Fig.1 of ref. [10]). The spin-flip kinetics within glauber kinetics algorithm defines the spin acceptance probability and also assumes ergodicity however the spin is chosen randomly [26,27]. The system quenches from a random disordered to ordered lattice phase and leads to formation of domains of various functional activities [7]. The Monte carlo simulations of the lattice for a particular coupling range n and synaptic strength J with T generated the temporal data for each spin (neuron). In order to substantiate the global system behavior we did average over temporal data of 10 ensemble sets in each case. The averaged time-series data is further used to construct binary undirected graphs by defining inter-neuronal connections using Pearson correlation coefficient for different coupling range and strengths (see Methods). Thus, our analysis of a 2D lattice of neurons with fixed position and randomly allocated spin or neuron activity at each site lead to the transformation into emergence of functional patterns of synchronous activity at the near-critical temperature as global synaptic strength.
Model based functional cortical patterns (FCPs). The emergence of spin clusters based on temporal correlations of neural activity, translated in the form of binary connectivity matrices, is the result of self-organized criticality that formulates the functional pattern of neurons in the brain [28]. These FCPs could mimic the functional networking that emerges in the brain cortex with taskspecific activation of neurons. The recurrence plots in Fig.1B shows the emergence of FCPs being more prominent at near critical temperature, T c for each coupling range. The long-ranged coupling at T = 1, is showing very less spin (neuron) clustering as the system is far from equilibrium and therefore it has many random connections not forming any ordered patterns or domains reflected in recurrence plots. However, at near-critical synaptic strength (T ∼ T c ), the transition to critical phase allows the formation of a number of functional patterns or domains. We could signify the near critical transition as a state of attaining the optimal synaptic strength to construct synchronous functional patterns of neural activity and that could mimic the conscious state of a simple cognitive task in the brain [6]. Further, at above-criticality (T > T c ) these patterns are violated due to enormous decline in randomness of neuronal connections, causing break-down of neuronal self-organization. Also, the high synaptic strength could signify the abnormal state with loss of functional connectivity as observed in diseases such as Alzheimer's etc [29]. The similar trend is followed by the critical and short coupling range (see in Fig.1B). However the underlying dynamics is different for short and long range interactions. The long-range coupling (n = 3) has been speculated with fast dynamical clustering at near critical temperature, T c = 2.9 J kB whereas, short-range coupling (n = 6) has exhibited slow dynamics at T c = 1.5 J kB [10]. From equation ( refTC), short range T c (n = 6) is given by , the long range T c (n = 3) is given by, This shows T c (N ) depends on system size, N , for long-range coupling only. These changes in T c could suggest the enlargement of FCPs, generated with respect to short and long range depending on N , with increasing synaptic strength signifying the local and global coupling over population of neurons [28,30].
Network topology and characterization. The complex network theory approach, originated from graph theory, delineates the consanguineous activity of neurons in terms of quantitative measures [31,32]. We have applied the network theory approach here to interpret the networks generated through the neuro-physiological EEG data set and the simulated dynamics of the population of neurons.

Properties of model based FCPs:
We have computed network theory attributes to study the topological properties of the FCPs generated (as shown in Fig.1), using our model approach. The sheer transition in the clustering coefficient distribution, C(k), plots from below (black) to near critical state (red) ( Fig. 2A column a) apparently classifies the organization of a functional cluster in a hierarchical pattern as exhibited through the power

FIG. 2:
Basic network attributes of model-based FCPs and EEG-based FNs: A) FCPs for the three distinguished coupling ranges at synaptic strengths below criticality (in black), near criticality (in red) and above criticality (in green). The columns exhibits a) clustering coefficient, C(k), versus degree (k) plots with degree exponent α. b) The probability of degree distribution, P(k), plots with degree exponent γ c) The Neighborhood connectivity, Cn(k), plots with degree exponent β. B) a) The connectivity trend of EEG data based FNs generated for three visual stimulus conditions S1 (single image shown), S2 match (two similar images) and S2 nomatch (two disimilar images) of one subject with 10 trials shown in different colors in each case. The columns b), c) and d) illustrates the C(k), P(k) and Cn(k) plots versus degree (k), respectively. The insider plots depict power law fitting trend of 10 subjects with their average degree exponents. Note -In all the coupling ranges long, critical and short their is transition shift at near criticality (shown in red) for all the three basic network attributes that characterizes the hierarchical and scale-free organization through the power law behaviour of their heavy-tailed distributions. The network measures are subsequently in similar trend for near-critical FCPs and EEG based FNs, apart from the obvious disparity among subjects.
law scaling behaviour C(k) ∼ k −α as a function of degree k with scaling exponent α for all the distinguished coupling ranges (long, critical and short). The data has been plotted on a logarithmic scale with non-linear curve fitting to expose the power-law scaling behavior. The cases with initial scattering of data points have shown that the power law is exhibited only by their long heavy tail whereas, other cases have shown power law scaling constituting whole or majority of data points. The power law distribution has been confirmed by the algorithm given by Clausset et al. [33]. However, at above critical state (green) due to increased in randomness of neuronal firing there are very less connections to form a functional pattern. Since α for long range coupling provides larger value (α = 0.58) as compared to short range coupling (α = 0.35), the long range coupling near critical temperature probably enhances the active self-organization of neurons and local effects are allowed to participate at global phenomena. The probability of degree distribution, P (k) again exhibits well defined fractal nature P (k) ∼ k −γ with scaling parameter γ near critical temperature with larger value of γ = 1.04 at long range coupling and shows similar transition into a power law function at near critical synaptic strength thus substantiates the presence of high degree nodes with less probability in the neural network. Similarly, the neighborhood connectivity distribution also follows power law nature, C n (k) ∼ k β with k, and near critical temperature exhibits positive β value for short range and critical couplings indicating possibility of formation of rich-club, where some of the hubs hand in hand together to control the network. However, long range coupling try to initiate neurons against this rich club formation (negative value in β). This assures that long range coupling tries to keep hierarchical organization in the functional cortical pattern near critical temperature phase. However, short and near critical range couplings try to bring the network towards near scale free features, where C n (k) becomes nearly independent of k and β becomes positive indicating significant importance of group of hubs (rich club formation) respectively. Further, the behaviour at below critical temperature exhibits positive β values (black colour data) for all types of coupling ranges (short, critical and long) which is the clear indication of rich club formation showing significant important role of hub groups in regulating FCPs. On the other hand, in the case of above criticality regime, proper and clear behaviour of network properties in the FCPs are not exhibited for all coupling ranges.

EEG-data based functional networks (FNs):
The EEG time series data of the whole brain while visualizing a particular single object (S1) , two matching objects (S2 match) and two different objects (S2 nomatch) in human subjects has been taken from the UCI database. Since the networks or visibility graphs constructed from the EEG time series data carry the properties embedded in the EEG data [12], characterization of these networks of various subjects may highlight how these brain functional networks work and organize themselves. Each node in the constructed networks corresponds to the time step mapped to the EEG time series data which could be the resultant of the algebraic sum of all possible interacting neuronal signals (including all possible coupling ranges) at that time step. Hence, these nodes in the networks could be thought of as the representations of the interacting neurons at various cortical regions of the brain at various time steps collected by the EEG probes as signals. Further, these functional networks (FNs) thus constructed from the EEG data of the subjects, while doing a visual task, exhibit well defined clusters of represented neurons with slight variation in their sizes and properties due to the different paradigms (S1, S2 match and S2 nomatch) ( Fig. 2B column a). These FNs of the subjects shown to specific visual stimulus task declares the dynamics of the neuronal functional patterns specific to the the visual cortex.
Now the connectivity plots and topological characterization of the EEG-based FNs show clear emergence of clusters/modules/domains of varied sizes in all subject cases (S1, S2 match and S2 nomatch) ( Fig.2B a). Similar to our theoretically obtained results ( Fig. 1), the behaviour of the topological properties of these networks (we present the results of only for 10 trials of single subject and its true for large number of trials of large number of subjects we have taken from the UCI database for analysis) follow power law nature indicating fractal nature of the system. More specifically, properties of our neuron activity pattern model near critical temperature is closely similar to the properties of FNs constructed from EEG data. Then, these topological parameters can be represented by, The positive value in exponent β in C n indicates the assortativity property in the network (Methods) which enable significant hubs in the network to coordinate among them [34]. Further, the emerged modules/clusters follow hierarchical features and hence have significant correlation among them (Methods). Because of these two reasons these FNs most likely exhibit rich-club formation which is consistent with the experimental report on brain [35]. This means that the rich-club formation of significant hubs in the hierarchically self-organized neuronal system may probably act to stabilize the system both at local and global levels.
Our interpretatory theory professes, as if a slow and weak signal (say at T < T c ) propagates among the cluster of neurons creating a pre-conscious state for neurons and then a perturbation (signal) would provide an optimal synaptic strength (T c )) that could drive the system towards a functional cognitive state. (This could give sense to the condition, when there are many processes or thoughts that passes through our mind in an unconscious state and an external stimulus is required to attain the conscious state.) The C(k) and P (k) distributions have been found to be in accordance with the brain functional connectivity network attributes based on the experimental data as reported in [6,36]. The results shown in [6] reflects the functional network attributes while listening to music and finger tapping. The complexity in functional brain networks has been classified as scale-free considering spatio-temporal neural activity through neuro-electrical and imaging techniques [37]. From this we could affirm the validity of our model approach in characterization of the emerging FCPs formed in the brain while some specific cognitive task is done [38]. Thus, the network representations of nodes at T ≃T c for long, critical and short range interactions assures the self-organization of nodes into a FCP which exhibits scale-invariant and hierarchical organization. This supports our understanding that the self-organization of the interacting neurons is maintained with the emergence of FCPs for disparate ranges at near criticality [30].
Centrality and rich-club measures: hierarchical to scale-free transition The roles and characterization of the significant hubs in the brain networks, in which emerged modules and sparsely distributed hubs work in self-organized manner, can be well studied and compare using centrality measurements and rich-club analysis of the networks constructed from both neuron pattern activity model and EEG data. The centrality measures, as shown in Fig.3, characterizes the FCPs, at below and near criticality for all coupling ranges and does a fair comparison with the FNs generated through the EEG data of a visual stimulus. Since the closeness centrality measurement C C of a node is estimated by the inverse of average distances of the node with the rest of the other remaining nodes [39], large value of C C , which is quantified by small average path distance, provides significant importance of the node in regulating the network. The C C as a function of k follows power law nature, C C (k) ∼ k δ , both in proposed model network and network constructed from EEG data except the value C C in networks of EEG data is comparatively smaller than that of model network (Fig. 3 (a)). The behaviour of C C for model network near critical temperature, T ∼ T c for critical interaction range is closely similar to that of network of EEG data, where, C C increases with k indicating significant role of high degree hubs in regulating brain network in terms of fast information processing of the neurons may be by coordinating the high degree nodes in all the modules of the networks. On the other hand, short and long range correlation of the neurons again trigger more importance of the hubs moving towards the scale-free phase where hubs are more important than modules in controlling the networks. For higher degree nodes, the value is less at T ≃T c , this confirms the rareness of hubs in the network.
Betweenness centrality of the networks, which is another centrality measure from the amount of paths passing through each node from the rest of the nodes in the network, quantify the volume of traffic diffusing through each node in these networks [39,40], and found that it follows power law behaviour, C B ∼ k λ in both model and EEG data networks (Fig. 3 (b)). Since C B increases as k increases, it clearly indicates that high degree nodes have high traffic of information, and they play crucial role to control and stabilize the network. In this case, the behaviour of the C B for all neuron coupling regimes (short, critical and long range) are closely mimic with that of networks of EEG data in qualitative sense (indicated by nearly parallel fitting lines on the data points). This nature of C B (k), characterizes the cluster at T ≃T c with higher value for high degree nodes and thus ensures greater significance of those nodes in fast information transmission across the network [32]. Similarly, the properties of the eigen-vector centrality, which is the measure of influence of a node to any other nodes in the network that could induce long-term traffic risk [41], show power law nature, C E (k) ∼ k ǫ to both for model and EEG data based networks indicating closely mimic behaviour (due to nearly parallel fitted lines). This eigenvector centrality, C E (k), further depicts the intensity of most prominent nodes in networks and found that network is highly communicative at T ≃T c with low cost. Hence, high degree nodes have significant influencing capability to other nodes as well as more risk of attack to them. These high degree nodes could be generally target nodes of any brain disorder (disease, functional disorder etc).
Now the behaviour of all three centrality measurements (C C , C B , C E ) indicate the significantly important roles of high degree nodes in both model and EEG data based networks in terms of fast signal processing, regulation of information traffic, influencing other nodes in the network. This could be a clear signature of rich-club formation in these networks, indicating significant role of the high degree hubs in regulating the networks but removal of these hubs will not cause network breakdown. Another thing depicts through our results is the presence of high degree nodes in long as well as short range connections. This is in accordance with the experimental results done by [42], where they found high degree neurons through direct recording from upto 500 neurons in sliced cultures of mouse somatosensory cortex. In all centrality measurements of model based FCPs at T ∼ T c and T < T c , long-range neuronal coupling contribution is higher than the remaining critical and short-range interactions indicating significantly more important role of long-range interaction of neurons in FCPs. One possible reason for more importance of long-range neuronal interaction in FCPs could be enhancement of propagation of localized properties/perturbation driven by short-range neuronal interaction at global level for better cross-talk to keep network stabilize in self-organized manner.
The plot in Fig. 3 d, quantifies the rich-club coefficient, R C , as a function of N k>k level , which is the number of nodes with degree (k) greater than a particular k level , where k level = 1, ..., max(k) (see Methods). This has encountered presence of hubs (high degree nodes) connecting module/communities patterned in rich-club arrangement, to maximize communication at low cost maintaining robustness of network. The comparative analysis has shown that FCPs at below critical strength are far from the behaviour of the empirical data networks used here. However, at near critical strength FCPs leads to a shift towards the behaviour of FNs generated. The degree exponents of various network attributes calculated has been shown in Table1 for all coupling range and strengths along with the EEG data. Apart from the visible differences in degree exponents of our model and EEG clusters, the trend followed by them in network attributes is quite similar.

Modular organization: Is near critical state active
Modularity signifies the formation of community or modular structure in the network [43]. The network centrality measures characterize the important nodes that determine the functional efficiency of the network. However, defining the role of a node based on its position within the module and participation among inter-linking modules is an effective way to portray complex networks [44]. In Fig.4A (a), we compare the community affiliated structure, (C i ), of long, critical and short-ranged coupling networks with the FNs of the EEG data. In column Fig. 4A (b), the within-module degree, Z i , versus participation coefficient (P i ) plots relates the organization of influential nodes based on their intra and inter-module links. The results show that the near-critical FCPs has marked good coherence with the EEG data. In Fig.4B, based on the categorization of nodes done in [44], we have shown the transition in role and position of module nodes in the near-critical modular networks. The P i versus degree k plots for the distinguished coupling ranges at below-critical and near-critical synaptic strengths has been shown in columns Fig.  4B (a) and (b) respectively. The nodes have been labelled as module hubs and non-hubs with Z i ≥ 2.5 and Z i < 2.5 respectively. The further segregation depending on their P i values describes nodes as, R1 ultra-peripheral non-hubs (P i ≤ 0.05) with all the edge connections in the same module; R2 peripheral non-hubs (0.05 < P i ≤ 0.62) with mostly intra-modular edges; R3 connector non-hubs (0.62 < P i ≤ 0.80) with many inter-modular edges; R4 kinless non-hubs (P i > 0.80) with homogeneous sharing of connections among modules; R5 provincial hubs (P i ≤ 0.30) with most of intra-modular connections; R6 connector hubs (0.30 < P i ≤ 0.75) with majority of inter-modular associations and R7 kinless hubs (P i > 0.75) with homogeneous associations among all the modules [44]. In our results of Fig.4B, the transition from temperature below to near-critical regime leads to the occurrence of non-hub connector nodes, R3, (shown in green) in all coupling ranges. The early emergence of such hub nodes in long range coupling could be due to its fast neuronal dynamics. Also, the near-critical FCPs has shown transition to optimum hierarchical organization of connector hub nodes, R6, and peripheral non-hub nodes, R2, with respect to degree for all coupling ranges, validating its functional modularity. The FNs of the EEG data (shown in violet), being a complex functional brain network, matches favorably with our model FCPs at near-criticality regime. Further in Fig. 5, we have analyzed the probability of degree distribution of nodes under a particular community, as constructed through the community affiliation Louvain method [45]. The modules generated through networks of below and near critical FCPs has been shown in Fig. 5 (a) and (b), in respective different color for each module. For comparative analysis the communities underlying the EEG data based networks has been shown in column (c) with quite similar trend among various trials (shown in respective different colors) of a single subject. The below to near critical transition signifies increase in the number of communities. Also, the near-critical communities at short-range organize themselves towards small-world and tend to become hierarchical and scale-free at long-range as depicted through their degree distributions. This could be a signature to differentiate local and global coupling mechanisms in the functional brain networks which sounds helpful to tackle specificity of brain disorders.

Multi-fractal signature due to complex brain organization
Fractal characterization of a complex temporal signal has been done using the multi-fractal detrended fluctuation analysis (MFDFA) [13,15]. The power law scaling behaviour exhibited by the networks at near-critical strength characterizes the approximately self-similar temporal pattern formation, which further announces the existence of a fractal dimension to describe the complex organization and role of short and long range correlation of the components of the system. The comprehensive methodology on M F DF A in [15,46] signifies the presence of multi-fractality in biological complex systems reflected in time-series data analysis of the system. We have studied the presence of multi-fractality in self-organizing FCPs at below and near critical strengths and compared the same with the EEG data generated FNs. In Fig.6 plot(a), the overall root mean square fluctuation function, F q , for q in range -5 to +5 is showing variation of fast and slow evolving fluctuations over segments of scale s in the time series, for each coupling range (long, short and critical) at T < T c and T ≃ T c and the EEG data as shown in respective colors. The scale of range 16 to 1024 has been used to calculate the function F q at each scale s and is same for both the model and empirical EEG data. The slope of this scaling function F q determines the q-order Hurst exponent H q . The H q is an important parameter to characterize multi-fractality as it captures the scaling behavior of small and large fluctuations in the time series data, thus, its dependence on q validates the multi-fractal nature. In Fig.6(b), the negative q exhibits the scaling trend of small fluctuations and positive q represents the scaling behaviour of large fluctuations in the data segments of particular scale. The Fig6(b) has shown variation mainly in the scaling nature of small fluctuations at negative q ′ s and depicts almost similar trend for large fluctuations at positive q ′ s. The case of EEG data demonstrates the result of one subject each of S1, S2 match and S2 nomatch whereas the model data is the averaged result of 10 trials of time-series for all the three cases (long, short and critical). The H q of small and large fluctuations in the model data might get affected because of the averaging of data. However, we are more concerned about the scaling trend of the near critical transition results being similar with the EEG data rather than the below criticality. The Fig6(c) represents the q-order mass exponent t q versus q, which exhibits the variation in mass exponent t q with respect to large (+ve q) and small (−ve q) fluctuations. The Fig6(d), partitioned into two subplots, further represents the generalized fractal dimension D q versus singularity exponent h q for below-criticality in the left subplot and near-criticality and EEG data in the right subplot. The mass exponent t q , is a factor used to define the fractal dimension. In Fig6(d), the clear transition in the h q and D q , from below to near-critical state, exhibits presence of multiple scaling exponents represented by the arc of multi-fractal spectrum of long, critical and short coupling ranges that signifies increase in complexity. The near-critical multi-fractal spectrum curves has shown There is emergence of R3 non-hub connector nodes and R7 kinless hubs in the near-critical FCPs. Also, the organization of R2 peripheral non-hub nodes and R6 connector hub nodes is showing the right hieararchy with respect to degree upon transition to near criticality. In a) and b) columns for long, critical and short range coupling we have shown probability of degree distribution, P km , of each node degree k in module m of the below and near critical FCPs, respectively. Here, different colors represent different communities. c) The column shows P km versus km plots of the communities underlying the EEG data based networks for the three paradigms S1, S2 match and S2 nomatch. Here, different colors represent communities of various trials of a single subject. Note-The modules of near-critical FCPs at short-range follows poisson distribution that transitions to the scale-free organization at long-range coupling.
high congruence with spectrum curves of the FNs of a visual task-based EEG data.

Conclusion
We have done the characterization of cortical neuronal circuitry in the form of FCPs, that emerges out in the brain during any cognitive assignment. We have analyzed combined effect of both coupling strength and coupling range as drivers defining the interactions in our model. Through graph theory measures, the effect of coupling range on the network of neural populations depicts clear transition in the dynamic state of connectivity at near-critical synaptic strength (T ≃ T c ). The network attributes calculated here, represents the topological information of the connected undirected networks generated for different coupling range. These results illustrate the transition of states from a randomly connected network featuring hierarchical properties to a state that owns scale free characteristics, exhibiting FCPs, with presence of hubs as shown by the heavy tail distribution of degree distribution plots. The probability of degree distribution P (k) and clustering coefficient C(k) plots clearly explain the emergent complex dynamics of functionally connected patterns (FCPs) that could mimic the functional neuronal patterns inside the brain. The dynamics of the networks at T ≃ T c for long, critical and short range couplings could correspond to the topology of emergent functional patterns formed in the brain while performing a specific cognitive task. The self-organization in the brain lead to the emergence of FCPs in the form of hierarchically organized functional modules as reflected from the power law nature of the rich-club organization and centrality measures of the networks constructed [38,47]. The presence of network centrality and rich-club behaviour, in similar trend with the EEG-based FNs, ensures the presence of coordinating hubs or dense modular cluster of neurons for effective communication in the model-based FCPs. This makes the functional pattern of neurons robust enough to maintain connectivity while wiring/rewiring among neurons as the dynamic pattern of connectivity among neurons affects the cost of information transmission [48]. The near-critical FCPs have shown presence of R6, connector hubs, and R7, kinless hubs (in Fig.4B), with inter-modular connections to broadcast the information. The FCPs of neurons exhibit the hierarchical and scale-free organization with distinct community/pattern formation, where, coordination of the hubs via cross-talking communities in the network is one of the key factors of information processing both at local and global levels of neuronal organization. In our study it is evident that this mechanism of FCPs organization is rectified by local and global perturbations in it, which can be propagated throughout the network with the help of various ranges of interaction (short, critical and long range interaction). This manifests the fractal properties of these FCPs, which further exhibited multifractal spectrum formation at near critical strength. Thus, the near-critical transition state marks the onset of cognitive state through self organizing dynamics at near criticality.
As we focussed on the basis of generation of these FCPs we made the following insights. Multiple stimulus attempts can cause increase in synaptic strength and make more long-range connections to the complex functional pattern in brain. The dynamics of these FCP's procreates diverse cognitive states in the brain. These cognitive states could be thought of as a bunch of meta-stable states in a system with near dynamic equilibrium. Brain dynamics is known to function at near criticality. From the study of near critical dynamics for long, critical and short-ranged functional patterns, we can interpret that criticality is a range and not just a point. This critical range could accommodate a set of dynamic functional states whose synchronous permutational amalgamation could lead its way to cognitive behaviors. Thus, self-organization is the key event that generates functional patterns in the brain cortex. Also, such analytical modeling of functional networks of neurons on the basis of inter-neuronal coupling range could further assist in extensive assessment of brain disorders, characterized by aberrant functional connectivity. The research in [17] have analyzed 2D ising model to exhibit topological characterization of functional networks of neuronal activity at criticality, being more robust against structural defects. However, in brain disorders, such as, ASD irregularity in the structural connectivity has caused less intricate long-ranged and more dense short-ranged functional connectivity [49]. Also, the extent in reflection of structural connectivity changes on the topology of functional network cannot be determined till date and sounds debatable. Though, the research in [17] makes our belief more strong towards the critical phase transition state as the emerging brain functional state. We hope our approach towards understanding the generation and characterization of functional patterns of neuronal activity in the brain cortex could embark us with new insights on brain functional irregularities.

Network Construction:
The functional connectivity of neuronal population in brain could be examined through a binary connectivity matrix where '1' represents connection between two nodes and '0' is antithesis of '1' [50]. We have computed the functional connectivity based on statistical consanguinity of neuronal time-series using Pearson-correlation coefficient, used to measure inter-neuronal correlations [50]. The Pearson correlation coefficient r xy , defines correlation between nodes x and y using the formula, Where, x i and y i represents temporal data of two spins x and y for the total number of monte carlo steps. The binary undirected matrix is transformed into graph G(S, E) consists of nodes set S = {s i | i ∈ 1, ..., N } where N being the total spin sites i.e. 4096 in our study. Then, based on the functional correlation of spins the delineated edge list is defined as E = {e ij | (s i , s j ) ∈ S × S}. We generated the binary connectivity matrix by averaging the data of 10 time-series from our model for each coupling range and global synaptic strength (Temperature). We applied a common thresholding limit for defining connections among the spins (neurons). We calculated the average correlation value in each case and then did the final average to get a unique threshold i.e. 0.2 for our model. This adjacency matrix is used to construct undirected networks for long, short and critical coupling range (n) at different overall synaptic strengths (T) employing NetworkX module in python. As we are studying the emergent behavior of our model generated functional patterns, we also constructed the functional brain cortical networks using the visibility graph approach on the EEG time-series data.
Visibility graph approach: The EEG data has been taken from the UCI database, where the subjects were shown to three types of visual stimulus S1 (single image shown), S2 match (two similar images) and S2 nomatch (two dissimilar images) and recording was done using 64 electrodes with sampling frequency of 256Hz for 1 second in healthy human subjects [14]. We applied the visibility graph approach in the EEG time-series data to construct the functional brain networks (FNs) generated while doing the visual task [12]. We took 10 trials for each type of stimulus. In this approach, each time state or neuronal state in the EEG time series (or neuron activity pattern model based time series) data, which is the resultant signal of the interacting neurons in the brain at that time state, is taken as node in the constructed network. The connections between ant two time or neuronal states with data value n a and n b at time point t a and t b respectively would be defined if the third time step n c (t c ) satisfies the following condition, The extracted network would be connected, at least to its neighbors, undirected and invariant according to the algorithm. The characteristic properties of the time-series get delineated in the form of resultant network.

Network theory attributes:
Network theory is the much applied theoretical approach by researchers for characterizing and studying complex brain networks [9]. We are providing some important attributes, which are used for our analysis, as given below.

Degree distribution
The probability of degree distribution to have degree k in a network defined by G (N, E), where N and E are sets of nodes and edges of the network, is given by, P (k) = n k N , where, n k and N are number of nodes having k degree and size of the network, respectively [51]. The degree distribution, P (k), for the Erdös-Reńyi random network follows Poisson distribution and deviates from the random for small-world networks, whereas, for scale-free and hierarchical networks it follows power law, P (k) ∼ k −γ , where, 2 < γ < 3 [52,53]. Clustering co-efficient Clustering co-efficient, which characterize how strongly a node is connected to the rest of the nodes in the network, of an ith node in the network can be estimated as the ratio of the number of its nearest neighborhood edges to the total possible number of edges of degree k i , C(k i ) = Ei k i C2 , where, E i and k i are the number of connected pairs of nearest-neighbors of ith node [51]. C(k) → constant for scale-free, random and small-world networks, whereas, for hierarchical network, C(k) ∼ k −α , with α ∼ 1 [52,53]. Neighborhood connectivity It is characterized by the measure of mean connectivities of the nearest neighbors of each node in a network given by, C n (k) = u uP (u|k), where, P (u|k) is the probability that a link belongs to a node with degree k points to a node with connectivity u [54,55]. Hierarchical network follows, C(k) ∼ k −β , with β ∼ 0.5 [55]. However, if C(k) ∼ k β (positive β), then the network exhibits assortative nature indicating possibility of coordinating high degree hubs in the network. Betweenness centrality Betweenness centrality of a node w is the measure of sharing amount that a node i needs w to reach j via shortest path, and is given by, C N (w) = i,j;i =j =w = dij (w) dij , where, d ij (w) is the number of geodesic paths from node i to j passing through w, and d ij indicates number of geodesic paths from node i to j [39]. It characterizes the amount of information traffic diffusing from each node to every other nodes in the network [40]. Closeness centrality Closeness centrality of a node u is characterized by the harmonic mean of the geodesic paths connecting u and any other nodes in the network, C C (u) = n i dui , where, d ui is the geodesic path length between nodes u and i, and n is total number of nodes connected to node u in the network [56]. It generally measures the rate of information flow in the network, where, larger the value of C C corresponds short path lengths, and hence fast information processing in the network, and vice versa [40]. Eigen-vector centrality Eigen-vector centrality of a node w in a network, which indicates the power of spreading of a particular node in the network, can be expressed by, C E (w) = 1 λ i=nn(w) u i , where, nn(w) is the nearest neighbors of node w, λ is the eigenvalue of eigenvector v i Av i = λv i , A is the adjacency matrix of of the network [56]. The principal eigenvector of A, which satisfy Au w = λ max , is taken to have positive eigenvector centrality score. Rich-club co-efficient Rich nodes corresponds to nodes having large links in the network, and they have tendency to form tight subgraph among them, which is referred to as rich-club formation in the network [57]. This rich-club phenomena can be quantified by rich-club coefficient defined by, R C (k) = E >k N >k C2 , where, E >k as the number of edges after removing nodes less than degree k distributed among N >k nodes [34]. It characterizes many important properties of the network, namely, information traffic backbone, mixing properties of the network, etc. The rich club organization analysis has been done to estimate presence of connector hubs. The Brain connectivity toolbox has been used to compute the rich club coefficient and network centrality measures [50,58].

Participation index
The participation coefficient (P i ) signifies the distribution of connections of a particular node i with respect to different communities [44] given by, where, k ic is the number of connections made by node i to nodes in module c and k i is the total degree of node i. The P i value determines the distributional uniformity of the neuronal connections, specified by the range 0-1. The escalating value signifies more homogeneous allocation of links among all the modules.
The within-module degree Z i is another measure to quantify the role of a particular node i in the module c i . High values of Z i indicate more intra-community connections than inter-community and vice-versa [44].
where, k i is the number of connections of the node i to other nodes in its module c i .k ci is the average k over all nodes in the module c i and σ kc i is the standard deviation of k in c i .

Community detection mehtod:
The complex network structure can be forbidden into communities or modules, specified with less than expected number of connections among them. We have used the community louvain method to design non-overlapping communities out of the network [45]. This method has outperformed other methods in terms of both modularity and computational time [45]. To create a significant division of a network the benefit function called modularity (Q) is defined as, The modularity Q, is maximized for good partitioning of the graph G(S, E) with N as total nodes. The A ij and B ij defines the exact and expected number of connections between nodes i and j. W = i,j A ij , c i and c j are the communities nodes i and j belongs to,δ ci,cj equals to 1, when nodes i and j falls in same community and 0, if they do not.

MFDFA analysis:
Multi-fractal detrended fluctuation analysis (MFDFA) approach has been used to ensure the fractality of dynamic patterns generated over non-stationary temporal data in complex biological systems [13,59]. The statistical fractals generated through physiological time series data show self affinity in terms of different scaling with respect to direction [59]. The presence of multiple scaling exponents can be examined in a time series using the matlab formulation of the MF-DFA method [46]. Parameters characterizing multifractality are scaling function (F), Hurst exponent (H), mass exponent (t), singularity exponent (h) and Dimension (D) as explained in [46]. For a time series signal x j of finite length l with random walk like structure, can be computed by the Root mean square (RMS) variation, where, x is the mean value of the signal, and i = 1, 2, ..., l. The signal X has been divided into n s = int( l s ) non-overlapping segments of equal size s. To avoid left-over short segments at the end, the counting has been done from both sides therefore 2n s segments are taken into account. This defines the scale (s) to estimate the local fluctuations in the time series. Thus, the overall RMS, F, for multiple scales can be computed using the equation, where, v = 1,2, ..., n s and x v (i) is the fitting trend in each segment v. The q-order RMS fluctuation function further determines the impact of scale (s) with large (+ q's) and small (-q's) fluctuations, as follows, The q-dependent fluctuation function F q (s) for each scale (s) will quantify the scaling behaviour of the fluctuation function for each q, F q (s) ∼ s Hq (15) where, H q is the generalized Hurst exponent, one of the parameter that characterizes multi-fractality through small and large fluctuations ( negative and positive q's) in the time series. The H q is related to q-order mass exponent t q as follows, from t q , the singularity exponent h q and Dimension D q is defined as, The plot of D q versus h q represents the multifractal spectrum of the time series.