# Observation of complex functional cortical patterns in brain cognition

Jasleen Gund1,2 and R.K. Brojen Singh1<EMAIL_ADDRESS>(Corresponding author)
1School of Computational and Integrative Sciences, Jawaharlal Nehru
University, New Delhi-110067, India.
2Cognitive Brain Dynamics Lab, National Brain Research Centre, Manesar-122050,
India.

###### Abstract

Synaptic plasticity and neuron cross-talk are some of the important key
mechanisms underlying formation of dynamic clusters of active neurons. The
essence of this study is to model and decipher the mechanism of emergence of a
task-specific functional connectivity among a population of neurons. We have
used our proposed neuron activity pattern (NAP) model to define an assorted
range of interactions and simulated across a 2D lattice with two-state
(firing/non-firing) neurons. We observed the emergence of scale-free,
hierarchically organized ordered patterns of active neurons, that we call as
functional cortical patterns (FCPs), only at the near-critical phase
transition. We have done extensive topological characterization of FCPs and
tested its congruency with the functional brain networks obtained from the
empirical electrophysiological data of a visual stimulus task. Our results of
network theory attributes, fractal analysis and functional cartography have
confirmed the transitioning effect while generation of FCPs. This gives us
strong implication that the functional neuronal system supports far-from-
equilibrium dynamics upon receiving a stimulus. Our investigation for the
long, critical and short-ranged interactions has lead to an interpretation
that near-critical transition phase is a spectrum of critical temperatures
(coupling strengths) being affected by the distance between the neurons. The
range of interaction among neurons plays an important role in defining the
functional severity of the neuronal circuitry. We propose that interaction
range as a function of coupling strength drives connectivity in an augmenting
fashion as we observed an increase in coupling strength from short to long-
range interactions, among the population of neurons. This gives us insight
towards the intensity of cognitive behaviour, which is attained by multiple
stimulus attempts.

Keywords: Long-range Ising model, Neuron activity patterns, task-specific
cognition, Brain functional networking, Near-critical dynamics.

##  I Introduction

The phenomenon of global features of complex cognitive functions in the brain
arise from the organization of basic functions systematically coordinated in
localized brain regions (Luria,, 1980). 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 (Chialvo,, 2010). 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.
Electroencephalogram (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 the concentration of neurotransmitter,
transfusion of ions through ion channels and voltage generated across the
membrane of pre and post-synaptic neurons (Klinshov et al.,, 2014; Voglis and
Tavernarakis,, 2006). These biological factors determine the strength of the
diffused signal among neurons and neuronal organization. Such recurring
signals strengthen 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 (Colicos
and Syed,, 2006; Eguíluz et al.,, 2005; Sporns et al.,, 2000, 2002; Sporns,,
2002). But, how the synchronous neuronal activation and generation of
functional patterns lead to the 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 how
this range is determined for information processing and shaping neural
circuitry among the neurons. 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 the system evolution McCoy and Wu, (2014). In our foregoing work, we
proposed the neuron activity pattern or NAP model which is based on the
kinetic Ising model system Gundh et al., (2015). In that work, we established
a critical interaction range that demarcates the long and short ranges of
interaction. Our approach evaluates the role of various interaction ranges on
the strength of coupling and the overall Hamiltonian dynamics, as compared to
the Onsager result that solved for only nearest-neighbour interactions
(Onsager,, 1944). It has been shown in (Gundh et al.,, 2015) that each
interaction/coupling range has their unique critical temperature for
generating spin clusters at the near-critical regime of phase transition
though, 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. However, the emergence of these functional
clusters and their organization are still an open question. We have termed
these patterns of neuronal activity as functional cortical patterns or FCPs.

In this study, we have generated model-based networks on the basis of,
firstly, different ranges of interaction (long, critical and short) and
secondly, another important parameter i.e. Temperature to define below, near
and above criticality regimes. We have taken this architecture, further, to
characterize temporal correlations based network connectivity of FCPs in our
model data and found its congruency with the task-specific brain functional
connectivity, attained through the empirical data. Since the functional
connectivity reflects statistical dependencies between spatially distributed
neuronal groups. We have made a comparison among networks of simulated and
empirical data based on their topological characterization. This study has
modelled the emergence of functional connectivity in the brain and to prove
its fidelity we have applied different approaches such as network theory
attributes, community detection, functional cartography of nodes and multi-
fractal analysis. We have used the EEG data of a specific visual task on
healthy human subjects for comparison with our model data (Begleiter,, 1999).
We have addressed important features of functional cluster formation, their
organization and properties at below, near and above criticality for each
interaction range. The major concern of this work is, whether the FCPs
generated through our neuron activity pattern model, could mimic the
topological characteristics of task-specific functional brain connectivity.

Figure 1: The Recurrence plots of model-based functional cortical patterns
(FCPs): A) Outline of the workflow used to characterize the model and
empirical time-series data and to attain congruency in their functional
activity patterns. B) For Long ($n=3$), Critical ($n=4$) and Short ($n=6$)
coupling ranges at synaptic strengths below-critical ($T=1$), near-critical
($T\sim T_{c}$) and above-critical ($T>T_{c}$) strengths ( as shown in black,
red and green respectively), in a network of 4096 nodes. Note:-  There is the
formation of communities of active neurons at near-criticality for all the
coupling range interactions.

##  II Results and Discussion

We first need to delineate the details of the model simulation in the next
sub-section, succeeded by the deliberation over results.

###  II.1 Correlated neuron activity could optimize brain function

The concept of self-organization rightly defends the dynamics of complex
systems (Kauffman,, 1993) and also aids in understanding brain functionality
(Bassett and Gazzaniga,, 2011). Thinking of neural dynamics from a statistical
mechanics window is a bit challenging but is ruling out since Ising model was
used to explain the dynamical avalanches of neuron activity in the brain
(Chialvo,, 2010; Eguíluz et al.,, 2005). 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 expound that the
brain works at a transition point when it receives the optimum signal strength
(Beggs and Timme,, 2012). Still, there is thirst for digging more into the
complex dynamics of brain functionality. How brain does the transition towards
a functional network that spans the cortical layers of brain, acquiring non-
random synaptic connectivity for performing a cognitive task (Song et al.,,
2005).

This research work has used the NAP model that was introduced in our previous
work (Gundh et al.,, 2015), therefore, the details regarding the Hamiltonian
dynamics of the Ising model and its numerical simulation using the Monte Carlo
Glauber kinetics with non-conserved magnetization has been briefly discussed
here. 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
(Fraiman et al.,, 2009; Goodarzinick et al.,, 2018). This makes us strongly
presume that phase transition state is the critical state of brain functional
dynamics at which cognitive actions take place. 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 contemplated the
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 result 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,

$\displaystyle H\>=-{\displaystyle{\displaystyle\sum_{<ij>}{\textstyle
J(r_{ij},n)s_{i}s_{j}},}}\;\quad\;\>\;s_{i}=\pm 1\,\forall\,i,$ (1)

Where $s_{i}$ and $s_{j}$ denotes the activity of neuron at site i and a
neighbouring site j and J represents the coupling strength among neurons
defined as a function of the distance between two $i$th and $j$th 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 behaviour given by,
$J(r_{ij},n)=\frac{J}{r_{ij}^{n}}$ (Cannas and Tamarit,, 1996), where, $d$ is
the dimension of the system, and $n$ can be expressed as $n=d+\sigma=2+\sigma$
(Fisher et al.,, 1972; Blanchard et al.,, 2013) which can qualify for short-
range for $\sigma>2$, whereas for long-range for $0<\sigma<2$ (Picco,, 2012;
Grüneberg and Hucht,, 2004). Further, $J>0$ always (ferromagnetic case) and
$n>0$ (Cannas and Tamarit,, 1996). We have taken a 2D square lattice of spins
with side $L$ and $N=L\times L$ being the total number of spin sites. To
define coupling range we have assumed the centre of the 2D square as the
reference 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 the thermodynamical parameters of the system,
one can define the following potential function $U(n)$ as a slowly decay
potential function (Chialvo,, 2010; Cannas and Tamarit,, 1996), where
$n=2+\sigma$ for the 2D system, scaling $J\rightarrow J/N$ in the limit
$N\rightarrow\infty$ and using Euler-McLaurin sum formula (Bruijn,, 1981), we
get the potential energy function for each interaction range n as follows,
(Gundh et al.,, 2015)

$\displaystyle U\left(n\right)$ $\displaystyle=$
$\displaystyle\underset{N\rightarrow\infty}{\lim}\frac{1}{N}\sum_{<i,j>}^{N}\frac{J}{r_{ij}^{\sigma}}\approx\lim_{N\rightarrow\infty}J\int_{1}^{\sqrt{N}}drg(r)r^{3-n}$
(2) $\displaystyle=$
$\displaystyle\lim_{N\rightarrow\infty}J\left[\begin{matrix}\frac{1}{2}\ln(N),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}n=4~{}(critical)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
\frac{1}{n-4}\left(1-N^{2-n/2}\right),~{}~{}for~{}~{}n>4~{}(short~{}range)~{}~{}~{}~{}~{}\\\
\frac{1}{4-n}N^{2-n/2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}0<n<4~{}(long~{}range)~{}\end{matrix}\right.$

where, g(r) is the pair distribution function such that $g(r)~{}1$ for $r>>1$.
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 (Hiley and Joyce,, 1965). Following this
technique one can estimate critical temperature from the numerically derived
closed form approximated equation (Hiley and Joyce,, 1965) given below,

$\displaystyle\frac{U(n)}{k_{B}T_{C}}=1+\frac{f_{2}}{U(n)^{2}}+O(U^{-4})$ (3)

where, $f_{2}=\sum_{<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\rightarrow
J/N$ in the limit $N\rightarrow\infty$. Now using Euler-McLaurin sum formula
(Bruijn,, 1981), $f_{2}$ is given by,

$\displaystyle f_{2}$ $\displaystyle=$
$\displaystyle\underset{N\rightarrow\infty}{\lim}\frac{1}{N}\sum_{<i,j>}^{N}\left[\frac{J}{r_{ij}^{\sigma}}\right]^{2}\approx\lim_{N\rightarrow\infty}J^{2}\int_{1}^{\sqrt{N}}drg(r)r^{5-2n}$
(4) $\displaystyle=$
$\displaystyle\lim_{N\rightarrow\infty}\frac{J^{2}}{2}\left[\begin{matrix}\left[1-\frac{1}{N}\right],~{}~{}~{}~{}~{}~{}for~{}~{}n=4~{}(critical)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
\frac{1-N^{3-n}}{n-3},~{}~{}~{}~{}~{}~{}for~{}~{}n>4~{}(short~{}range)~{}~{}~{}~{}~{}\\\
\frac{N^{3-n}-1}{3-n}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}0<n<4~{}(long~{}range)~{}~{}\end{matrix}\right.$

Now using equations (3), (2) and (4) we could able to arrive at an expression
for finite critical temperature given below,

$\displaystyle T_{c}(n,N)$ $\displaystyle=$
$\displaystyle\frac{1}{k_{B}}\frac{1}{1+\frac{f_{2}}{U(n)^{2}}}$ (5)
$\displaystyle=$
$\displaystyle\frac{J}{k_{B}}\left[\begin{matrix}\lim_{N\rightarrow\infty}\frac{N[\ln(N)]^{3}}{2N+N[\ln(N)]^{2}-1},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}n=4~{}(Critical)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
\frac{4(n-3)}{(n-2)^{2}(n-4)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}n>4~{}(short~{}range)~{}~{}~{}~{}~{}\\\
\lim_{N\rightarrow\infty}\frac{2}{4-n}\frac{N^{3(4-n)/2}}{2N^{4-n}+(4-n)^{2}\left[\frac{N^{3-n}-1}{3-n}\right]}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}for~{}~{}0<n<4~{}(long~{}range)~{}~{}\end{matrix}\right.$

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 a
modified phase. On the other hand, for fixed $n$, $T_{c}$ increases as $N$
increases both in critical and long-range regime. The quantities such as range
of interaction $n$ (or coupling range) and coupling strength $J$ controls the
order parameter i.e. temperature $T$ of a spin-lattice system. 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})\approx 1$ (Onsager,, 1944), 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}\rightarrow[2.269-2.29]$ (Binder and Heermann,, 1988). The reason
for variation in the critical parameter value could be due to finite size
($N$) of the system (Binder and Heermann,, 1988), 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. (Gundh et
al.,, 2015)), signified by the sudden drop in the total magnetization (order
parameter) of the system.

Our idea states that this modelling framework mimics the underlying dynamics
among a population of neurons. This study intends to relate the terminology
being used in two different yet relatable systems, to characterize the brain
as a dynamic physical system. We have defined the temperature, analogously, as
the global synaptic strength given to a system of neurons (or groups of
neurons) that generates the synchronous simultaneous activation of neurons
sharing common synaptic strength. Therefore, the critical temperature would
mirror the optimal perturbation (stimulus) imparted to the system i.e. relayed
and associated in the form of synaptic strength in the brain. We have taken a
2D $64\times 64$ lattice system of 4096 spin sites depicting either firing or
non-firing neuron state, allocated randomly, and simulated for 500 monte carlo
steps at a defined global synaptic strength. 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\frac{J}{k_{B}}$ as below critical temperature for all coupling ranges.
However, the near-critical transition temperature, $T_{c}$, is different for
the long $(T_{c}\sim 2.9\frac{J}{k_{B}})$, critical $(T_{c}\sim
1.9\frac{J}{k_{B}})$ and short $(T_{c}\sim 1.4\frac{J}{k_{B}})$ ranged
interactions, based on their respective dynamics (see Fig.1 of ref. (Gundh et
al.,, 2015)). The spin-flip kinetics within Glauber kinetics algorithm defines
the spin acceptance probability and also assumes ergodicity, however, the spin
is chosen randomly (Binder and Heermann,, 1988; Newman and Barkema,, 1999).
The system quenches from a random disordered to ordered lattice phase and lead
to the formation of domains/patterns of functional activity (Sporns et al.,,
2000). The Monte Carlo simulations for a particular coupling range $n$ and
global synaptic strength $T$ has generated the temporal activity data based on
the state of each neuron. In order to substantiate the global system
behaviour, 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 the 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 activity at each site lead to the emergence of FCPs of synchronous
activity at the near-critical temperature, the optimum synaptic strength.

Figure 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 exhibit a) clustering coefficient, C(k), versus degree (k)
plots with degree exponent $\alpha$. b) The probability of degree
distribution, P(k), plots with degree exponent $\gamma$ c) The Neighborhood
connectivity, $C_{n}(k)$, plots with degree exponent $\beta$. 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 dissimilar images) of one subject with 10 trials shown in
different colours in each case. The columns b), c) and d) illustrates the
C(k), P(k) and $C_{n}(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 characterize 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.

###  II.2 Model-based functional cortical patterns

The neuron activity patterns can be imprinted in a connectivity matrix based
on the correlation in the activity of neurons. The temporal correlations of
the neural activity has been translated in the form of binary connectivity
matrices. We have applied different approaches for building 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 have tagged the model-based functional
clusters of synchronous activity as FCPs. On the other hand, the visual task-
based EEG time-series has been transformed into a functional brain network,
FN, using the visibility graph approach (Lacasa et al.,, 2008). We have
acquired significant congruence in the FCPs and FNs in their topological
properties and functional characterization by applying the network theory and
multi-fractal approaches (Sporns,, 2002; Kantelhardt et al.,, 2002). The model
simulations at various interaction ranges and temperature strengths
substantiate the emergence of functional patterns of activity only at the
near-critical strength as depicted through the recurrence plots in Fig.1B. In
this work, we have extended the interpretation of our model data and applied a
fixed global synaptic strength which determines the intensity of overall
coupling among the system of N neurons, being classified as below-criticality,
near-criticality and above-criticality of the phase transition dynamics. The
critical temperature $T_{c}$, which signifies the state of near-critical phase
transition, is found to mimic with the optimum synaptic strength to coordinate
short and long-range neurons in the brain for the generation of functional
clusters of neurons at local and global levels.

The emergence of FCPs in our model is the result of self-organized criticality
that formulates the functional pattern of neurons in the brain (Kitzbichler et
al.,, 2009). These FCPs could mimic the functional networking that emerges in
the brain cortex with task-specific 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 clustering and therefore it has many
random connections not forming any ordered patterns or domains as reflected in
recurrence plots. However, at near-critical synaptic strength ($T\sim 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 (Eguíluz et al.,, 2005).
Further, at above-criticality ($T>T_{c}$) these patterns are violated due to
enormous decline in the 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 (Stam et al.,, 2007). 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\frac{J}{k_{B}}$ whereas,
short-range coupling ($n=6$) has exhibited slow dynamics at
$T_{c}=1.5\frac{J}{k_{B}}$ (Gundh et al.,, 2015). From equation (5), short-
range $T_{c}$ $(n=6)$ is given by $T_{c}=\frac{3}{8}\frac{J}{k_{B}}$. Whereas,
since $\lim_{n\rightarrow 3}\frac{N^{3-n}-1}{3-n}=\lim_{n\rightarrow
3}\frac{-(\ln(N))N^{3-n}}{-1}=\ln(N)$, the long-range $T_{c}$ $(n=3)$ is given
by,

$\displaystyle T_{c}(N)\approx\frac{J}{k_{B}}\frac{2N^{3/2}}{2N+\ln(N)}$ (6)

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 (Kitzbichler et al.,, 2009; Deco et al.,, 2012).

###  II.3 Network topology and characterization.

The complex network theory approach, originated from graph theory, delineates
the consanguineous activity of neurons in terms of quantitative measures (Papo
et al.,, 2014; Bullmore and Sporns,, 2009). We have applied the network theory
attributes here to characterize the functional networks generated through the
neuro-physiological EEG data set and the simulated FCPs of the population of
neurons.

####  II.3.1 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
network attributes calculated here, represents the topological information of
the connected undirected networks generated for different coupling ranges of
interaction. 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-law scaling behaviour
$C(k)\sim k^{-\alpha}$ as a function of degree $k$ with scaling exponent
$\alpha$ 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 behaviour. The cases with an 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. (Clauset et al.,,
2009). However, at the above critical state (green) due to increased in the
randomness of neuronal firing there are very fewer connections to form a
functional pattern. Since $\alpha$ for long-range coupling provides larger
value ($\alpha=0.58$) as compared to short-range coupling ($\alpha=0.35$), the
long-range coupling at near-critical temperature probably enhances the active
self-organization of neurons and local effects are allowed to participate in
global phenomena. The probability of degree distribution, $P(k)$ again
exhibits well defined fractal nature $P(k)\sim k^{-\gamma}$ at near-critical
temperature with larger value of scaling parameter, $\gamma=1.04$, for long-
range coupling and shows similar transition into a power-law function, thus,
substantiate the presence of high degree nodes having less occurrence
probability in the neural network. Similarly, the neighbourhood connectivity
distribution also follows power-law nature, $C_{n}(k)\sim k^{\beta}$ with $k$,
and near-critical temperature exhibits positive $\beta$ 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 $\beta$). 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 $\beta$ becomes positive indicating
significant importance of a group of hubs (rich club formation) respectively.
Further, the behaviour at below critical temperature exhibits positive $\beta$
values (black colour data) for all types of coupling ranges (short, critical
and long) which is a clear indication of rich club formation showing the
significantly important role of hub groups in regulating FCPs. On the other
hand, in the case of the above criticality regime, proper and clear behaviour
of network properties in the FCPs are not exhibited for all coupling ranges.

####  II.3.2 Properties of EEG based functional networks:

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 (Lacasa et al.,, 2008),
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 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 (in Fig. 1), the behaviour of the topological
properties of these networks (we presented the results only for 10 trials of a
single subject and its also true for all trials of the subjects we have taken
from the UCI database for analysis) follow power-law nature indicating the
fractal nature of the system. More specifically, the 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,

$\displaystyle\bf{\Gamma(k)}\sim
k^{F};~{}~{}\bf{\Gamma(k)}=\left[\begin{matrix}P(k)\\\ C(k)\\\
C_{n}(k)\end{matrix}\right],~{}~{}F=\left[\begin{matrix}\gamma\\\ \alpha\\\
\beta\end{matrix}\right]\rightarrow\left[\begin{matrix}-2.16-2.31\\\
-0.73-0.8\\\ 0.21-0.23\end{matrix}\right]$ (7)

The positive value in exponent $\beta$ in $C_{n}$ indicates the assortativity
property in the network (Methods) which enable significant hubs in the network
to coordinate among them (Colizza et al.,, 2006). Further, the emerged
modules/clusters follow hierarchical features and hence have a 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 the brain (Teller et al.,, 2014). 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 interpretation 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 (Eguíluz et al.,, 2005; Sporns et al.,, 2004). The results shown
in (Eguíluz et al.,, 2005) 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 (Beggs and Plenz,,
2003). 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(Thomas Yeo et al.,, 2015). Thus, the network
representations of nodes at $T\,\simeq$$T_{c}$ for long, critical and short-
range interactions assures the self-organization of nodes into an 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
(Deco et al.,, 2012).

Figure 3: Comparative analysis of network attributes: The model-based FCPs
have been analysed to mark the near-critical phase transition in comparison
with the EEG data based FNs. a) The closeness centrality, $C_{C}(k)$, versus
degree (k) plot with degree exponent $\delta$. b) The betweenness centrality,
$C_{B}(k),$ plot with degree exponent $\lambda$. (C) The eigen-vector
centrality, $C_{E}(k)$, plot with degree exponent $\epsilon$. (D) The rich-
club coefficient, $R_{C}(k)$, is plotted against $N_{k>k_{level}}$i.e. number
of nodes with degree (k) greater than particular $k_{level}$. Note\- The near-
critical FCPs have shown the right transition shift towards the FNs of a
visual task-based data. For a particular degree node the $C_{B}$ and $C_{E}$
value has increased to make its role more economical, though not much change
in $C_{C}$values observed. The $R_{C}(k)$ plot has shown power-law trend with
heavy tail distribution.

###  II.4 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 a 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 (Linton C. Freeman,, 1978), a 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)\sim k^{\delta}$, 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 the
model network (Fig. 3 (a)). The behaviour of $C_{C}$ for model network at
near-critical temperature, $T\sim T_{c}$ for critical interaction range is
closely similar to that of the network from EEG data, where, $C_{C}$ increases
with $k$ indicating a 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\,\simeq$$T_{c}$, this confirms the rareness of hubs in
the network.

Betweenness centrality of the networks, which is another centrality measure
from the number 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 (Linton C. Freeman,, 1978; Borgatti,, 2005), and found that it
follows power-law behaviour, $C_{B}\sim k^{\lambda}$ in both model and EEG
data networks (Fig. 3 (b)). Since $C_{B}$ increases as $k$ increases, it
indicates that high degree nodes have high traffic of information, and they
play a 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\,\simeq$$T_{c}$ with
higher value for high degree nodes and thus ensures greater significance of
those nodes in fast information transmission across the network (Bullmore and
Sporns,, 2009). 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 (Bonacich,, 1972), show power-law
nature, $C_{E}(k)\sim k^{\epsilon}$ 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\,\simeq$$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}$)
indicates 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 the
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 are the presence of high degree nodes in long as well as
short-range connections. This is in accordance with the experimental results
done by (Nigam et al.,, 2016), where they found high degree neurons through
direct recording from up to 500 neurons in slice cultures of mouse
somatosensory cortex. In all centrality measurements of model based FCPs at
$T\sim T_{c}$ and $T<T_{c}$, long-range neuronal coupling contribution is
higher than the remaining critical and short-range interactions indicating a
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 enhanced propagation of localized
properties/perturbation, driven by short-range neuronal interaction, at the
global level for better cross-talk to keep network stabilization in a self-
organized manner.

The plot in Fig. 3d, 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 a rich-club arrangement, to maximize
communication at low-cost maintaining robustness of the 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 have
been shown in Table-1 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 are quite
similar.

###  II.5 Modular organization: Is near-critical state active

Modularity signifies the formation of community or modular structure in the
network (Newman,, 2006). 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 (Guimera and Amaral,, 2005). 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 relate 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 (Guimera and Amaral,, 2005), 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}\geq 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}\leq 0.05$) with all the edge connections in the
same module; $R2$ peripheral non-hubs ($0.05<P_{i}\leq 0.62$) with mostly
intra-modular edges; $R3$ connector non-hubs ($0.62<P_{i}\leq 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}\leq 0.30$)
with most of the intra-modular connections; $R6$ connector hubs
($0.30<P_{i}\leq 0.75$) with majority of inter-modular associations and $R7$
kinless hubs ($P_{i}>0.75$) with homogeneous associations among all the
modules (Guimera and Amaral,, 2005). In our results of Fig.4B, the transition
of temperature from 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
the optimum hierarchical organization of connector hub nodes, $R6$, and
peripheral non-hub nodes, $R2$, with respect to the 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 favourably with
our model FCPs at the near-criticality regime. Further in Fig. 5, we have
analyzed the probability of the degree distribution of nodes under a
particular community, as constructed through the community affiliation Louvain
method (Blondel et al.,, 2008). 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 have been shown in column
(c) with quite a similar trend among various trials (shown in respective
different colors) of a single subject. The below to near-critical transition
signifies an 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.

Figure 4: Community structure in networks of the model and EEG data: A) The
model-based FCPs for the long, critical and short-ranged coupling has been
compared with the EEG data based FNs (shown in violet). a) The community
affiliation index, ($C_{i}$), versus degree (k) plots. b) The within-module
degree, $Z_{i}$, versus participation coefficient, $P_{i}$, plots exhibits
presence of module-hubs based on the classification scheme defined as R1-
ultra peripheral non-hubs, R2- peripheral non-hubs, R3- connector non-hubs,
R4- kinless non-hubs, R5- provincial hubs, R6- connector hubs, R7- kinless
hubs. B) The classification of module nodes, based on the $Z_{i}$ and $P_{i}$
values, and their organization depicted through the $P_{i}$ versus degree (k)
plots in the model generated FCPs for a) Below critical synaptic strength at T
= 1.0 and b) Near critical strength of their respective $T_{C}$. The presence
of various hub and non-hub nodes has been coded in respective colors shown in
the figure. Note- The congruency of the EEG data (shown in violet) with the
near-critical dynamics (shown in red) claims the model robustness. 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 hierarchy with respect to degree upon
transition to near criticality. Figure 5: Analysis of Probability of degree-
distributions in modules of networks: In a) and b) columns for long, critical
and short-range coupling we have shown probability of degree distribution,
$P_{k_{m}}$, 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_{k_{m}}$ versus $k_{m}$ 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.

###  II.6 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) (Kantelhardt et al.,,
2002; Song et al.,, 2005). 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 $MFDFA$ in (Song et al.,, 2005; Ihlen,, 2012) 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\simeq 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 a particular scale.
The Fig6(b) has shown variation mainly in the scaling nature of small
fluctuations at negative $q^{\prime}s$ and depicts an almost similar trend for
large fluctuations at positive $q^{\prime}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 the multi-fractal spectrum of long, critical and short coupling
ranges that signifies an increase in complexity. The near-critical multi-
fractal spectrum curves have shown high congruence with spectrum curves of the
FNs of a visual task-based EEG data.

Figure 6: Multi-fractal detrended fluctuation analysis (MFDFA) : The technique has been applied to the model-based FCPs at below and near-critical synaptic strengths and compared with the EEG-based FNs to determine multi-fractal nature. a) The q-order RMS fluctuation function $F_{q}$versus scale $s,$ quantifies the scaling trend in fluctuations. b) The q-order Hurst exponent $H_{q}$, determines the scale of the fluctuation function $F_{q}$. c) The q-order mass exponent $t_{q}$ d) The q-order singularity exponent $h_{q}$ versus Dimension $D_{q}$ exhibits the multifractal spectrum. Note- The transition of FCPs at near-critical state is clear in the $D_{q}$ versus $h_{q}$ plot that exhibits a perfect arc of multifractal spectrum around the range of exponents and complements with the FNs of a visual task-based EEG data. Degree exponents: | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $\epsilon$ | $\lambda$ | $\kappa$
---|---|---|---|---|---|---|---
FCPs ($T\langle T_{c}$) | Long | -0.32 | -0.03 | 0.1 | 0.1 | 1.04 | 1.99 | -0.15
Critical | -0.15 | -0.01 | 0.42 | 0.2 | 0.99 | 1.98 | -0.08
Short | -0.34 | -0.07 | 0.14 | 0.26 | 0.96 | 2.0 | -0.10
FCPs (T $\simeq$ $T_{c}$) | Long | -0.58 | -0.34 | -1.04 | 0.0 | 0.98 | 2.16 | -1.11
Critical | -0.4 | -0.07 | -0.6 | 0.08 | 0.88 | 1.97 | -0.72
Short | -0.35 | 0.0 | -0.6 | 0.01 | 0.98 | 1.93 | -0.78
EEG FNs | S1 | -0.8 | 0.21 | -2.31 | 0.10 | 1.6 | 2.92 | -0.7
S2 match | -0.73 | 0.23 | -2.1 | 0.10 | 1.8 | 2.34 | -0.78
S2 nomatch | -0.73 | 0.21 | -2.1 | 0.11 | 1.56 | 2.15 | -0.87
Table 1: The measure of degree exponents for various network theory attributes
(as mentioned in section 4.3).

##  III Conclusion

We have made an attempt to characterize cortical neuronal circuitry that
emerges out in the brain during any cognitive assignment. In our model, the
two-state dynamics has shown ordered pattern dynamics with the emergence of
functional clusters, when it evolves at near-critical regime for all the
interaction ranges. The emergent network connectivity at near-critical
strength, that we call as FCPs, symbolizes the emergence of functional
connectivity in the brain. We have done extensive topological characterization
of the FCPs and compared with the task-specific functional brain network of
human subjects. In this comprehensive work, we have validated our insights of
brain functionality by characterizing properties of self-organized functional
connectivity embedded in the time-series data of task-specific activity. 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 population depicts clear
transition in the dynamic state of connectivity at the near-critical strength
($T\simeq T_{c}$). This illustrates the transition from a randomly connected
network to a hierarchical network that owns scale-free characteristics with
the 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, explains the emergent complex
dynamics of FCPs. The dynamics of the networks at $T\simeq T_{c}$ for long,
critical and short-range coupling is found congruent with the emergent
functional connectivity 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 (Thomas Yeo et al.,, 2015; Meunier et al.,, 2010). 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 (De Domenico et
al.,, 2016). The near-critical FCPs have shown the 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.

Brain dynamics is known to function at near-criticality. There are recent
studies that depicts the non-linearity of the functional brain dynamics at the
near-critical regime (Tagliazucchi et al.,, 2016; Cocchi et al.,, 2017; Ezaki
et al.,, 2020). Our results of network theory attributes, fractal analysis and
functional cartography have confirmed the generation of a scale-free
hierarchical network (FCPs) only at the near-critical phase transition. We
have made the following insights focussing on the basis of generation of these
FCPs. The dynamics of these FCP’s could procreate 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. The near-critical
regime could accommodate a set of dynamic functional states whose synchronous
permutational amalgamation could lead its way to cognitive behaviours. Also,
from our study of near-critical dynamics for long, critical and short-ranged
functional patterns, we have interpreted that criticality is a range and not
just a point. As we found that the long-ranged dynamics require high coupling
strength to attain criticality than the short-ranged. This explains how the
increase in synaptic strength due to multiple stimulus attempts would make
more long-range connections to the complex functional pattern in the brain.
The dynamicity of long-range connections defines the task-specific functional
connectivity at the global level (Park and Friston,, 2013). Also, such
analytical modelling of functional networks of neurons based on inter-neuronal
coupling range could further assist in the extensive assessment of brain
disorders, characterized by aberrant functional connectivity. The research in
(Goodarzinick et al.,, 2018) 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 denser short-ranged functional connectivity
(Conti et al.,, 2017). 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
(Goodarzinick et al.,, 2018) makes our belief more strong towards the critical
phase transition state as the emerging functional state in the brain. We
believe that our approach towards understanding the generation and
characterization of functional patterns of neuronal activity in the brain
cortex embarks us with new understanding on brain functional irregularities
and also motivates modelling of functional connectivity in brain disorders.

There are obvious limitations that our model dynamics is binary and exhibits
only firing/non-firing states of a neuron or group of neurons at a single
lattice site. It overlooks the biochemical fluxes of a nerve cell and rather
focuses on the global connectivity state. In our study, we have considered
$64\times 64$ lattice system for our model and the empirical data is from 64
electrodes over the brain, this surely sounds a huge gap. However, brain
activity is known to follow statistical fractal nature, which allows the
characteristic properties to remain similar over the size (Franca et al.,,
2018). Also, we have addressed analytically that with an increase in system
size the critical temperature will increase in the long-range regime, but this
would not affect its near-critical dynamics. In spite of the huge complexity
gap, we have got significant congruency in terms of the characteristic power
law trends of the network topology between the task-specific functional brain
networks and the FCPs at the near-critical regime only. This research has
given us strong implication that the functional neuronal system supports far-
from-equilibrium dynamics upon receiving a stimulus that leads to a second-
order phase transition and self-organization is the key event that generates
functional patterns in the brain cortex. This work gives us very useful
insights towards the mechanistic view of the emerging brain functional
connectivity, though there is still much more to dig and modelling through
such simplistic model aid in understanding the concepts from a coarse lens
that surely can be investigated to microscopic levels.

##  IV Methodology

###  IV.1 Network Construction:

The functional connectivity of neuronal population in the brain could be
examined through a binary connectivity matrix where ’1’ represents the
connection between two nodes and ’0’ is antithesis of ’1’ (Sporns,, 2007). We
have computed the functional connectivity based on statistical consanguinity
of neuronal time-series using Pearson-correlation coefficient, used to measure
inter-neuronal correlations (Sporns,, 2007). The Pearson correlation
coefficient $r_{xy}$, defines correlation between nodes $x$ and $y$ using the
formula,

$\displaystyle
r_{xy}=\frac{\sum(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum(x_{i}-\bar{x})^{2}}\sqrt{\sum(y_{i}-\bar{y})^{2}}}$
(8)

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}\mid i\>\in
1$, …, $N$} where $N$ is 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}\mid(s_{i},s_{j})\in S\times 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.

###  IV.2 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 the sampling frequency of 256Hz for 1 second in
healthy human subjects (Beggs and Timme,, 2012). 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 (Lacasa et al.,, 2008).
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,

$\displaystyle\frac{n_{b}-n_{c}}{t_{b}-t_{c}}>\frac{n_{b}-n_{a}}{t_{b}-t_{a}}$
(9)

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.

###  IV.3 Network theory attributes:

Network theory is the much applied theoretical approach by researchers for
characterizing and studying complex brain networks (Sporns,, 2002). 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)=\frac{n_{k}}{N}$, where, $n_{k}$ and $N$ are number of nodes
having $k$ degree and size of the network, respectively (Newman,, 2009). The
degree distribution, $P(k)$, for the Erd$\ddot{o}$s-Re$\acute{n}$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)\sim k^{-\gamma}$, where, $2<\gamma<3$ (Albert and Barabási,,
2002; Barabasi and Oltvai,, 2004).
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})=\frac{E_{i}}{{}^{k_{i}}C_{2}}$, where, $E_{i}$ and $k_{i}$ are the
number of connected pairs of nearest-neighbors of ith node (Newman,, 2009).
$C(k)\rightarrow constant$ for scale-free, random and small-world networks,
whereas, for hierarchical network, $C(k)\sim k^{-\alpha}$, with $\alpha\sim 1$
(Albert and Barabási,, 2002; Barabasi and Oltvai,, 2004).
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)=\sum_{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$ (Maslov and Sneppen,, 2002; Pastor-
Satorras et al.,, 2001). The hierarchical network follows, $C_{n}(k)\sim
k^{-\beta}$, with $\beta\sim 0.5$ (Pastor-Satorras et al.,, 2001). However, if
$C_{n}(k)\sim k^{\beta}$ (positive $\beta$), then the network exhibits
assortative nature indicating the 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_{B}(w)=\sum_{i,j;i\neq j\neq w}=\frac{d_{ij}(w)}{d_{ij}}$, 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$ (Linton C. Freeman,, 1978). It characterizes the amount of information
traffic diffusing from each node to every other node in the network
(Borgatti,, 2005).
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)=\frac{n}{\sum_{i}d_{ui}}$, 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 (Canright and Engø-Monsen,, 2004). 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 (Borgatti,, 2005).
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)=\frac{1}{\lambda}\sum_{i=nn(w)}u_{i}$, where, $nn(w)$ is the nearest
neighbors of node $w$, $\lambda$ is the eigenvalue of eigenvector $v_{i}$
$Av_{i}=\lambda v_{i}$, $A$ is the adjacency matrix of the network (Canright
and Engø-Monsen,, 2004). The principal eigenvector of A, which satisfy
$Au_{w}=\lambda_{max}$, is taken to have positive eigenvector centrality
score.
Rich-club co-efficient
Rich nodes correspond to nodes having large links in the network, and they
tend to form tight subgraph among them, which is referred to as rich-club
formation in the network (Zhou and Mondragón,, 2004). This rich-club
phenomenon can be quantified by rich-club coefficient defined by,
$R_{C}(k)=\frac{E_{>k}}{{}^{N_{>k}}C_{2}}$, where, $E_{>k}$ as the number of
edges after removing nodes less than degree $k$ distributed among $N_{>k}$
nodes (Colizza et al.,, 2006). 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
the presence of connector hubs. The Brain connectivity toolbox has been used
to compute the rich-club coefficient and network centrality measures (Sporns,,
2007; Rubinov and Sporns,, 2010).
Participation index
The participation coefficient ($P_{i}$) signifies the distribution of
connections of a particular node i with respect to different communities
(Guimera and Amaral,, 2005) given by,

$\displaystyle P_{i}=1-\sum_{c=1}^{N}\left(\frac{k_{ic}}{k_{i}}\right)^{2},$
(10)

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 a 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 (Guimera and
Amaral,, 2005).

$\displaystyle Z_{i}=\frac{k_{i}-\bar{k_{c_{i}}}}{\sigma_{k_{c_{i}}}},$ (11)

where, $k_{i}$ is the number of connections of the node i to other nodes in
its module $c_{i}$. $\bar{k_{c_{i}}}$ is the average k over all nodes in the
module $c_{i}$ and $\sigma_{k_{c_{i}}}$ is the standard deviation of k in
$c_{i}$.

Community detection method:
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 (Blondel et al.,, 2008). This method has outperformed other
methods in terms of both modularity and computational time (Blondel et al.,,
2008). To create a significant division of a network the benefit function
called modularity (Q) is defined as,

$\displaystyle
Q=\frac{1}{W}\sum_{i,j=1}^{N}[A_{ij}-B_{ij}]\delta_{c_{i},c_{j}},$ (12)

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=\underset{i,j}{\sum}A_{ij}$, $c_{i}$and $c_{j}$are the communities nodes i
and j belongs to,$\delta_{c_{i},c_{j}}$equals to 1, when nodes i and j falls
in same community and 0, if they do not.

###  IV.4 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 (Kantelhardt et al.,, 2002; Eke et
al.,, 2002). The statistical fractals generated through physiological time-
series data show self affinity in terms of different scaling with respect to
direction (Eke et al.,, 2002). The presence of multiple scaling exponents can
be examined in a time-series using the Matlab formulation of the MF-DFA method
(Ihlen,, 2012). Parameters characterizing multifractality are scaling function
(F), Hurst exponent (H), mass exponent (t), singularity exponent (h) and
Dimension (D) as explained in (Ihlen,, 2012). 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, $X_{i}=\sum_{j=1}^{i}(x_{j}-\left\langle
x\right\rangle)$, where, $\left\langle x\right\rangle$ is the mean value of
the signal, and $i=1,2,...,l$. The signal X has been divided into
$n_{s}=int(\frac{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,

$\displaystyle
F^{2}(s,v)=\frac{1}{s}\sum_{i=1}^{s}\left\\{X[(v-1)s+i]-x_{v}(i)\right\\}^{2}$
(13)

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,

$\displaystyle
F_{q}(s)=\left\\{\frac{1}{2n_{s}}\sum_{\nu=1}^{2n_{s}}[F^{2}(v,s)]^{\frac{q}{2}}\right\\}^{\frac{1}{q}}$
(14)

The q-dependent fluctuation function $F_{q}(s)$ for each scale (s) will
quantify the scaling behaviour of the fluctuation function for each q,

$\displaystyle F_{q}(s)\thicksim s^{H_{q}}$ (15)

where, $H_{q}$ is the generalized Hurst exponent, one of the parameters that
characterize 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,

$\displaystyle t_{q}=qH_{q}-1$ (16)

from $t_{q}$ , the singularity exponent $h_{q}$ and Dimension $D_{q}$ is
defined as,

$\displaystyle h_{q}=\frac{dt_{q}}{dq},\quad\quad D_{q}=qh_{q}-t_{q}$ (17)

The plot of $D_{q}$ versus $h_{q}$ represents the multifractal spectrum of the
time-series.

##  V Conflict of Interest Statement

The authors declare no competing financial interests.

##  VI Author Contributions

R.K.B.S. and J.G. conceived the model. J.G. did the numerical experiment and
prepared the figures of the numerical results. J.G. and R.K.B.S. analyzed and
interpreted the results, and wrote the manuscript. All author read and
approved the manuscript.

##  VII Funding

J.G. and R.K.B.S. are financially supported by Council of Scientific and
Industrial Research through sanctioned project 25(0221)/13/EMR-II.

##  VIII Acknowledgments

We acknowledge the funding agencies and Jawaharlal Nehru University. This
manuscript has been released as a Pre-Print at
[https://www.biorxiv.org/content/10.1101/569244v1] (Gund and Singh,, 2019).

##  IX Data Availability Statement

The dataset analyzed for this study can be found in the [UCI machine learning
repository: EEG database data set
[https://archive.ics.uci.edu/ml/datasets/EEG+Database].

## References

  * Albert and Barabási, (2002) Albert, R. and Barabási, A.-L. (2002).  Statistical mechanics of complex networks.  Reviews of modern physics, 74(1):47.
  * Barabasi and Oltvai, (2004) Barabasi, A.-L. and Oltvai, Z. N. (2004).  Network biology: understanding the cell’s functional organization.  Nature reviews genetics, 5(2):101.
  * Bassett and Gazzaniga, (2011) Bassett, D. S., and Gazzaniga, M. S. (2011).  Understanding complexity in the human brain.  Trends in cognitive sciences, 15(5):200–209.
  * Beggs and Plenz, (2003) Beggs, J. M. and Plenz, D. (2003).  Neuronal Avalanches in Neocortical Circuits.  The Journal of Neuroscience, 23(35):11167–11177.
  * Beggs and Timme, (2012) Beggs, J. M. and Timme, N. (2012).  Being Critical of Criticality in the Brain.  Frontiers in Physiology, 3(June):1–14.
  * Begleiter, (1999) Begleiter, H. (1999).  UCI machine learning repository: EEG database data set.
  * Binder and Heermann, (1988) Binder, K. and Heermann, D. W. (1988).  Monte Carlo Simulation in Statistical Physics.  Springer Berlin Heidelberg.
  * Blanchard et al., (2013) Blanchard, T., Picco, M., and Rajabpour, M. a. (2013).  Influence of long-range interactions on the critical behavior of the Ising model.  EPL (Europhysics Letters), 101(5):56003.
  * Blondel et al., (2008) Blondel, V., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008).  Fast unfolding of communities in large networks.  Journal of Statistical Mechanics: Theory and Experiment, 2008(10):P10008.
  * Bonacich, (1972) Bonacich, P. (1972).  Factoring and weighting approaches to status scores and clique identification.  Journal of mathematical sociology, 2(1):113–120.
  * Borgatti, (2005) Borgatti, S. P. (2005).  Centrality and network flow.  Social networks, 27(1):55–71.
  * Bruijn, (1981) Bruijn, N. d. (1981).  Asymptotic methods in analysis.  Dover Publications, Inc. New York.
  * Bullmore and Sporns, (2009) Bullmore, E. and Sporns, O. (2009).  Complex brain networks: graph theoretical analysis of structural and functional systems.  Nature reviews. Neuroscience, 10(3):186–198.
  * Cannas and Tamarit, (1996) Cannas, S. A. and Tamarit, F. A. (1996).  Long-range interactions and nonextensivity in ferromagnetic spin models.  Physical Review B, 54(18):R12661–R12664.
  * Canright and Engø-Monsen, (2004) Canright, G. and Engø-Monsen, K. (2004).  Roles in networks.  Science of Computer Programming, 53(2 SPEC. ISS.):195–214.
  * Chialvo, (2010) Chialvo, D. R. (2010).  Emergent complex neural dynamics.  Nature Physics, 6(10):744–750.
  * Clauset et al., (2009) Clauset, A., Shalizi, C. R., and Newman, M. E. (2009).  Power-law distributions in empirical data.  SIAM review, 51(4):661–703.
  * Cocchi et al., (2017) Cocchi, L., Gollo, L. L., Zalesky, A., and Breakspear, M. (2017).  Criticality in the Brain: A synthesis of neurobiology, models and cognition.  Progress in neurobiology, 158:132–152.
  * Colicos and Syed, (2006) Colicos, M. a. and Syed, N. I. (2006).  Neuronal networks and synaptic plasticity: understanding complex system dynamics by interfacing neurons with silicon technologies.  The Journal of experimental biology, 209(Pt 12):2312–2319.
  * Colizza et al., (2006) Colizza, V., Flammini, A., Serrano, M. A., and Vespignani, A. (2006).  Detecting rich-club ordering in complex networks.  Nature, 2(February):110–115.
  * Conti et al., (2017) Conti, E., Mitra, J., Calderoni, S., Pannek, K., Shen, K. K., Pagnozzi, A., Rose, S., Mazzotti, S., Scelfo, D., Tosetti, M., Muratori, F., Cioni, G., and Guzzetta, A. (2017).  Network over-connectivity differentiates autism spectrum disorder from other developmental disorders in toddlers: A diffusion MRI study.  Human Brain Mapping, 38(5):2333–2344.
  * De Domenico et al., (2016) De Domenico, M., Sasai, S., and Arenas, A. (2016).  Mapping multiplex hubs in human functional brain networks.  Frontiers in Neuroscience, 10(JUL):1–14.
  * Deco et al., (2012) Deco, G., Senden, M., and Jirsa, V. (2012).  How anatomy shapes dynamics: a semi-analytical study of the brain at rest by a simple spin model.  Frontiers in Computational Neuroscience, 6(September):1–7.
  * Eguíluz et al., (2005) Eguíluz, V. M., Chialvo, D. R., Cecchi, G. A., Baliki, M., and Apkarian, A. V. (2005).  Scale-free brain functional networks.  Physical Review Letters, 94(1):1–4.
  * Eke et al., (2002) Eke, a., Herman, P., Kocsis, L., and Kozak, L. R. (2002).  Fractal characterization of complexity in temporal physiological signals.  Physiological measurement, 23(1):R1–R38.
  * Ezaki et al., (2020) Ezaki, T., dos Reis, E. F., Watanabe, T., Sakaki, M., and Masuda, N. (2020).  Closer to critical resting-state neural dynamics in individuals with higher fluid intelligence.  Communications Biology, 3(1):1–9.
  * Fisher et al., (1972) Fisher, M. E., Ma, S.-k., and Nickel, B. G. (1972).  Critical Exponents for Long-Range Interactions.  Physical Review Letters, 29(14):917–920.
  * Fraiman et al., (2009) Fraiman, D., Balenzuela, P., Foss, J., and Chialvo, D. R. (2009).  Ising-like dynamics in large-scale functional brain networks.  Phys. Rev. E, 79(6):61922.
  * Franca et al., (2018) Franca, L. G. S., Miranda, J. G. V., Leite, M., Sharma, N. K., Walker, M. C., Lemieux, L., and Wang, Y. (2018).  Fractal and multifractal properties of electrographic recordings of human brain activity: toward its use as a signal feature for machine learning in clinical applications.  Frontiers in physiology, 9:1767.
  * Goodarzinick et al., (2018) Goodarzinick, A., Niry, M. D., Valizadeh, A., and Perc, M. (2018).  Robustness of functional networks at criticality against structural defects.  Physical Review E, 98(2):022312.
  * Grüneberg and Hucht, (2004) Grüneberg, D. and Hucht, A. (2004).  Universal finite-size scaling analysis of Ising models with long-range interactions at the upper critical dimensionality: Isotropic case.  Physical Review E, 69(3):036104.
  * Guimera and Amaral, (2005) Guimera, R. and Amaral, L. (2005).  Functional catography of complex metabolic networks.  Nature, 433:895–900.
  * Gundh et al., (2015) Gundh, J., Singh, A., and Singh, R. K. B. (2015).  Ordering Dynamics in Neuron Activity Pattern Model: An Insight to Brain Functionality.  PLoS.ONE., 10(1932-6203 (Electronic)):e0141463.
  * Gund and Singh, (2019) Gund, J. and Singh, R. K. B. (2019).  Emergence of Functional Cortical Patterns of neurons characterize the self-organizing way to cognition in brain.  https://www.biorxiv.org/content/10.1101/569244v1
  * Hiley and Joyce, (1965) Hiley, B. J. and Joyce, G. S. (1965).  The Ising model with long-range interactions.  Proceedings of the Physical Society, 85(3):493.
  * Ihlen, (2012) Ihlen, E. A. F. (2012).  Introduction to multifractal detrended fluctuation analysis in Matlab.  Frontiers in Physiology, 3 JUN(June):1–18.
  * Kantelhardt et al., (2002) Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., and Stanley, H. E. (2002).  Multifractal detrended fluctuation analysis of nonstationary time series.  Physica A, 316(1-4):87–114.
  * Kauffman, (1993) Kauffman, S. A. (1993).  The origins of order: Self-organization and selection in evolution.  OUP USA.
  * Kitzbichler et al., (2009) Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009).  Broadband Criticality of Human Brain Network Synchronization.  PLoS Computational Biology, 5(3):e1000314.
  * Klinshov et al., (2014) Klinshov, V. V., Teramae, J.-n., Nekorkin, V. I., and Fukai, T. (2014).  Dense Neuron Clustering Explains Connectivity Statistics in Cortical Microcircuits.  PLoS ONE, 9(4):e94292.
  * Lacasa et al., (2008) Lacasa, L., Luque, B., Ballesteros, F., Luque, J., and Nuno, J. C. (2008).  From time-series to complex networks: The visibility graph.  Proceedings of the National Academy of Sciences, 105(13):4972–4975.
  * Linton C. Freeman, (1978) Linton C. Freeman (1978).  Centrality in Social Networks Conceptual Clarification.  Social Networks, 179(1968):215–239.
  * Luria, (1980) Luria, A. R. (1980).  Higher cortical functions in man.  Basic Books.
  * Maslov and Sneppen, (2002) Maslov, S. and Sneppen, K. (2002).  Specificity and stability in topology of protein networks.  Science, 296(5569):910–913.
  * McCoy and Wu, (2014) McCoy, B. M., and Wu, T. T. (2014).  The two-dimensional Ising model.  Dover publications: New York, USA.
  * Meunier et al., (2010) Meunier, D., Lambiotte, R., and Bullmore, E. T. (2010).  Modular and hierarchically modular organization of brain networks.  Frontiers in Neuroscience, 4(DEC):1–11.
  * Newman, (2009) Newman, M. (2009).  Networks:An Introduction.  OUP Oxford.
  * Newman and Barkema, (1999) Newman, M. and Barkema, G. (1999).  Monte carlo methods in statistical physics chapter 1-4.  Oxford University Press: New York, USA.
  * Newman, (2006) Newman, M. E. J. (2006).  Modularity and community structure in networks.  Proceedings of the National Academy of Sciences of the United States of America, 103(23):8577–8582.
  * Nigam et al., (2016) Nigam, S., Shimono, M., Ito, S., Yeh, F.-C., Timme, N., Myroshnychenko, M., Lapish, C. C., Tosi, Z., Hottowy, P., Smith, W. C., Masmanidis, S. C., Litke, A. M., Sporns, O., and Beggs, J. M. (2016).  Rich-Club Organization in Effective Connectivity among Cortical Neurons.  Journal of Neuroscience, 36(3):670–684.
  * Onsager, (1944) Onsager, L. (1944).  Crystal statistics. i. a two-dimensional model with an order-disorder transition.  Physical Review, 65(3-4):117.
  * Papo et al., (2014) Papo, D., Buldu, J. M., Boccaletti, S., and Bullmore, E. T. (2014).  Complex network theory and the brain.  Philosophical Transactions of the Royal Society B: Biological Sciences, 369(1653):20130520–20130520.
  * Park and Friston, (2013) Park, H. J., and Friston, K. (2013).  Structural and Functional Brain Networks: From Connections to Cognition.  Science, 342(6158):1238411–1238411.
  * Pastor-Satorras et al., (2001) Pastor-Satorras, R., Vázquez, A., and Vespignani, A. (2001).  Dynamical and correlation properties of the internet.  Physical review letters, 87(25):258701.
  * Picco, (2012) Picco, M. (2012).  Critical behavior of the Ising model with long-range interactions.  arXiv preprint arXiv:1207.1018.
  * Rubinov and Sporns, (2010) Rubinov, M. and Sporns, O. (2010).  Complex network measures of brain connectivity: Uses and interpretations.  NeuroImage, 52(3):1059–1069.
  * Song et al., (2005) Song, S., Sjöström, P. J., Reigl, M., Nelson, S., and Chklovskii, D. B. (2005).  Highly nonrandom features of synaptic connectivity in local cortical circuits.  PLoS biology, 3(3):e68.
  * Sporns, (2002) Sporns, O. (2002).  Network analysis, complexity, and brain function.  Complexity, 8(1):56–60.
  * Sporns, (2007) Sporns, O. (2007).  Brain connectivity.  Scholarpedia, 2(10):4695.
  * Sporns et al., (2004) Sporns, O., Chialvo, D. R., Kaiser, M., and Hilgetag, C. C. (2004).  Organization, development and function of complex brain networks.  Trends in Cognitive Sciences, 8(9):418–425.
  * Sporns et al., (2000) Sporns, O., Tononi, G., and Edelman, G. M. (2000).  Theoretical Neuroanatomy: Relating Anatomical and Functional Connectivity in Graphs and Cortical Connection Matrices.  Cerebral Cortex, 10(2):127–141.
  * Sporns et al., (2002) Sporns, O., Tononi, G., and Edelman, G. M. (2002).  Theoretical neuroanatomy and the connectivity of the cerebral cortex.  Behavioural Brain Research, 135(1-2):69–74.
  * Stam et al., (2007) Stam, C. J., Jones, B. F., Nolte, G., Breakspear, M., and Scheltens, P. (2007).  Small-world networks and functional connectivity in Alzheimer’s disease.  Cerebral Cortex, 17(1):92–99.
  * Tagliazucchi et al., (2016) Tagliazucchi, E., Chialvo, D. R., Siniatchkin, M., Amico, E., Brichant, J. F., Bonhomme, V., and Laureys, S. (2016).  Large-scale signatures of unconsciousness are consistent with a departure from critical dynamics.  Journal of The Royal Society Interface, 13(114):20151027.
  * Teller et al., (2014) Teller, S., Granell, C., De Domenico, M., Soriano, J., Gómez, S., and Arenas, A. (2014).  Emergence of assortative mixing between clusters of cultured neurons.  PLoS computational biology, 10(9):e1003796.
  * Thomas Yeo et al., (2015) Thomas Yeo, B. T., Krienen, F. M., Eickhoff, S. B., Yaakub, S. N., Fox, P. T., Buckner, R. L., Asplund, C. L., and Chee, M. W. L. (2015).  Functional specialization and flexibility in human association cortex.  Cerebral Cortex, 25(10):3654–3672.
  * Voglis and Tavernarakis, (2006) Voglis, G. and Tavernarakis, N. (2006).  The role of synaptic ion channels in synaptic plasticity.  EMBO reports, 7(11):1104–10.
  * Zhou and Mondragón, (2004) Zhou, S. and Mondragón, R. J. (2004).  The rich-club phenomenon in the internet topology.  IEEE Communications Letters, 8(3):180–182.