{'fraction_non_alphanumeric': 0.028555813635742586, 'fraction_numerical': 0.011989342806394316, 'mean_word_length': 4.67704091526081, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 11, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling brain dynamics by Boolean networks

Understanding the relationship between brain architecture and brain function is a central issue in neuroscience. We modeled realistic spatio-temporal patterns of brain activity on a human connectome with a Boolean networks model with the aim of computationally replicating certain cognitive functions as they emerge from the standardization of many fMRI studies, identified as patterns of human brain activity. Results from the analysis of simulation data, carried out for different parameters and initial conditions identified many possible paths in the space of parameters of these network models, with normal (ordered asymptotically constant patterns), chaotic (oscillating or disordered) but also highly organized configurations, with countless spatial–temporal patterns. We interpreted these results as routes to chaos, permanence of the systems in regimes of complexity, and ordered stationary behavior, associating these dynamics to cognitive processes. The most important result of this work is the study of emergent neural circuits, i.e., configurations of areas that synchronize over time, both locally and globally, determining the emergence of computational analogues of cognitive processes, which may or may not be similar to the functioning of biological brain. Furthermore, results put in evidence the creation of how the brain creates structures of remote communication. These structures have hierarchical organization, where each level allows for the emergence of brain organizations which behave at the next superior level. Taken together these results allow the interplay of dynamical and topological roots of the multifaceted brain dynamics to be understood.

Modelling rationale. The need to have a detailed experimental characterization of the interactional kinetics of the nodes at the brain connectome level, carefully studying the integration of both the structural and the functional connectome, at different scales of the neural organization, considering all the above-mentioned issues is gaining momentum [38][39][40][41][42] .
But the question of how connectivity between brain areas emerges remains underexplored. To fill this gap, a new modelling of brain networks, by using established connectomes and Boolean Networks (BN) is presented. By using the connectome as a structural network, and applying the BN approach to model the interactional dynamics of the nodes, we analyze the resulting behavior as functional neuronal dynamics. Brain dynamics can be described by many mathematical models, both discrete (such as Cellular Automata [43][44][45][46][47] , Boolean networks 48 , Probabilistic Boolean networks 49 ) that continuous [such as Cellular Neural Networks (CNNs) 50 , Ordinary Differential Equations (ODEs) or Partial Differential Equations (PDEs) 51,52 ]. We opted against using PDEs in multiple dimensions because they involve local dynamics, which are ill-suited for modeling the above-described connectome networks. By using sets of one-dimensional PDEs also introduces difficulties. While such systems efficiently describe the transmission of signals on a line, they do not provide specific information on how to explain what happens at the intersection of multiple lines, which occur as singular points. The introduction of Cellular Neural Networks (CNN) systems, by making time continuous but discretizing space, overcomes this problem, reducing network links to transmission lines [53][54][55][56] . Such modeling enables not only modeling the functioning of the brain areas of interest, which support various cognitive functions, but also to provide information about the ability to convey information across transmission lines. But, as we are interested in this work in how the process of generating cognitive circuits emerges to study the brain areas involved that give rise to specific cognitive functions, the best choice seemed that of adopting an extremely simple conceptualization such as the Boolean Networks models, which could form the baseline for further and more complex studies. Furthermore, by using BNs versus Cellular automata, one can indeed consider influences on the cognitive dynamics that are not determined by the immediate neighbors. It is obvious that thinking of the state of individual regions with only 2 values instead of a continuous state as in CNNs, or considering the system deterministic in its interactions instead of probabilistic, severely limits the possible dynamics. Despite this, the behaviors exhibited by the system, as we shall see later, maintain their complexity, typical of computational systems that, beyond the simplicity of the model, exhibit complex or chaotic behaviors and universal computational capabilities 57 . Thus, the Boolean networks model seemed to us to be the most suitable for analyzing the dynamic behavior of networks of this type without losing complexity and without introducing additional arbitrary elements into the modeling.
There are almost 6 different models of Boolean networks, which can be suitably exploited to achieve different dynamics, within the context of brain networks. a. CRBNs (Classic Random Boolean Networks), characterized by a synchronous pattern and a deterministic updating: in the passage from time t to time t + 1 all nodes are updated simultaneously, according to the Boolean functions associated with them. These rules are not totalistic as each node evolves according to the Boolean function that characterizes it. b. ARBNs (Asynchronous Random Boolean Networks) is an asynchronous not deterministic method. At each step, a node, and only one, randomly selected, is chosen and updated with probability 1/N. c. DARBNs (Deterministic Asynchronous Random Boolean Networks) is an asynchronous and deterministic method where two parameters are assigned to each node p and q (after choosing a maximum value for p and q ), with q < p . Each node will be updated in the r th step if q = r(modp). If more nodes follow the updating www.nature.com/scientificreports/ condition in one step of evolution, they are updated one by one, taking into account the changes already made during this step. d. GARBNs (Generalized Asynchronous Random Boolean Networks) is a semi-synchronous and non-deterministic method. After choosing a g number of nodes, at each step only those are updated. In this case, each node will be updated with probability equal to g. e. DGARBNs (Generalized Deterministic Asynchronous Random Boolean Networks) is a semi-synchronous and deterministic method. It follows the same rule DARBN, but if more nodes comply with the updating condition, they are updated simultaneously. f. PBN (Probabilistic Boolean Networks) 49 .
Consequently, in this paper we propose to model neural dynamics, maintained across a variety of different brain configuration states, at different levels of temporal and spatial scales, considering different sets of structural connectomes, from normal to many neurological disorders functioning 58 .
By using the Boolean networks approach, varying the initial conditions of the connectome pattern, we can explore how the functional interplay among brain regions operate at the macroscopic scale, with the aim to discover how the information flows in the underlying structural network, and how the dynamics arises from networks of neurons at the microscopic scale. From the simulation runs, plentiful rich dynamics emerged, giving much information on the above-mentioned issues and prospecting interesting results.
The paper is organized as follows: in Section "Methods", we provide a review of the methods used in this work, from the approach to map on the human connectome cognitive functions related to relevant fMRI studies, to that of the implemented Boolean Network model, to the analysis of emergent circuits as large-scale structures; Section "Results" deals with the results obtained for the adopted methods; Section "Circuit analysis" considers the neuronal circuits emerging from the simulations and their correspondence to the chosen cognitive reference model; Section "Influence of the connection threshold on the stability of the emerging circuits" explores the influence of threshold parameters on the stability of emerging circuits; finally, a discussion on practical implications drawn from this approach for the network neuroscience closes the work.

Methods
To model cognition with Boolean Networks, the methodological pipeline represented in Fig. 1 was adopted.
Step 1: mapping cognition into the human connectome. In this paper, we simulated the brain activation points for the Action-Execution behavioral domain as identified by Smith et al. 59 and stored in the Brain-Map.org site. They were mapped on the 82 ROI connectome network (shaped on the Brodmann Atlas) and on the network developed by Mingrui et al. 60 .
Basically, we started by downloading the activation maps for the cognitive function of Action-Execution, by the behavioral taxonomy of Brainmap.org 59,61,62 , now the principal database of fMRI and PET brain activation studies. Smith et al. 41 , identified the major functional networks of the brain, representative of a significant portion of all functional activation studies performed on more than 7000 functional maps, derived from different experimental fMRI conditions, summarized in terms of coordinated positions of activation peaks, using independent component analysis (ICA), a powerful data-driven approach for finding independent patterns in multivariate data stored in the BrainMap.org data repository. We chose to simulate Action-Execution behavior, a very complex activity involving a very large neural network. After extracting the brain activation points for the Action-Execution activation map, we mapped these points with corresponding networks nodes in the brain connectome. To obtain the brain connectome 18 , we used the classic methods related to the acquisition of DTI and high-resolution T1-weighted MRI, the Segmentation of the brain's white matter, the creation of whitematter tractography and the consequent segmentation into regions of interest (ROI), followed as standard in many software applications 63 , while the degree of an ROI allowed the definition of the role that each module plays as a source of connections with other ROI, thus improving the modelling process and the attainment of computational dynamics by using the Boolean networks 58 . In Network Neuroscience research 64 , a connectome is a mathematical representation of the brain by using nodes and edge, where each node is a region of the brain and the edges the functional links connecting the brain. The node degree corresponds to the number of edges that are attached to each node, that is the brain connectivity. Networks can be decomposed into communities or modules. Connections (edges) are either linking nodes within modules or between modules. Highly connected nodes are hubs, and they either connect primarily with other nodes in the same community (provincial hub) or with nodes that belong to different communities (connector hub). The mapping of physical and functional connections between subdivided anatomical areas of the brain with the ICAs categories of behavioral domain is a challenging problem since there are different brain parcellization methods which are related to cyto-architecture such as Brodmann 65 and others drawn from the Talairach Atlas 66 , based on different numbers of ROIs which have been embodied in software tools 60 . Though Brodmann's reference system of cortical maps does not match contemporary anatomical and functional data in many brain regions 67 and new approaches are required, it is still used to register functional activations of anatomical structures 68 . For this reason, we used an 82-connectome model as the parcellization method, based on Brodmann areas, mapping the ICAs of the behavioral domain directly on the ROIs of the Brodmann organization. This allowed us to identify the network node corresponding to the functional activation areas for each category of the chosen behavioral domain, obtaining the reference areas in the 82-connectome model, involved in the process under analysis. In Table 1, Brodmann areas with the corresponding nodes of the connectome has been represented within the Tailarach space for the mapping of each node into the ICAs behavioral domain area. www.nature.com/scientificreports/ Areas have been connected by an adjacency matrix of 82 rows and 82 columns, whose values are comprised between 0 and 1, with an average connection value of 0.51 and a standard deviation of 0.288 for each node. This means that, on average, each network node is strongly connected with the other nodes. From a modeling point of view, since it is known that if the nodes in a network are highly connected, the system does not exhibit emergent and/or complex behaviors, going directly onto fixed-point behavior or dying immediately, we chose to set the nodes' weight values at 0.87.
Step 2: Boolean network modelling. The concerns previously expressed about the necessity to make use of models that are explicitly defined in a computational framework is strictly tied to the requirements of accurately choosing the variables to be the most congruent possible with brain data and to explicitly and quantitatively parametrize them explicitly and quantitatively 61 . For satisfying this requirement, we have chosen a Boolean Network (BN). In the form of random Boolean Networks (RBN), the system was introduced long ago 69 , with the aim to study the regulation of gene networks. In this pioneering method, each gene had two states "0" and "1", where the value "1" stands for the gene activation state, on the contrary the value "0" stands for an inhibition state. Since then, many studies have used this model 70 , in the sector of systems biology 71 , especially for its ability to intercept the gene networks dynamics. The status of the nodes in a Boolean network changes each time step in relation to the status of the nodes connected to it and to the Boolean function that is assigned. The number of possible states of a Boolean network, with N nodes, is equal to 2 N . This number represents the cardinality of the state space. Boolean networks can be considered synchronous or asynchronous, deterministic or nondeterministic, based on their updating rule table. In synchronous dynamics, all nodes update their status at the same time creating a "massively parallel process". In cases where the state of the system can be predictable with certainty, the dynamics of the system is deterministic. The development rules can also be totalistic in the sense that they might be the same for all nodes of the network; otherwise, they can vary according to the node. These features are very useful for sophisticated modeling of the brain cognitive functions. Formally, a Boolean network is a triplet {G, S, f } , composed of a graph G , a set S = {0, 1} and an evolution rule f : C → C , where C is the space of all the possible configurations. As is well known, G = {V , E} is made by the set V = {v 1 , v 2 , · · · .v n } of n elements called vertices and by the set E = {ε 1 , ε 2 , . . . , ε s } of any edges between two vertices ε i ≡ e jk ≡ v j → v k . www.nature.com/scientificreports/ Each vertex v i can take a Boolean value of state s i = 0 or s i = 1 . The critical setting of the initial state will produce different emerging configurations of the system. This is particularly suitable for simulating brain dynamics. A configuration C = {s 1 , s 2 , . . . s n } is given by the values of the states of all the nodes. The cardinality of the set C is equal to 2 n . Given a node v i , the set of all nodes v k for which there is an ε m ∈ E such that e ki = v k → v i or e ik = v i → v k , it is said neighbor of v i , and it is indicated with I i . If there is an e m = v k → v i , then there is also an e n = v i → v k . In this condition, we have an undirected graph or just a graph. The set of all nodes v k is called the input neighborhood and is indicated with I i i , while the set of all nodes v k is named output neighbor and is indicated with , the graph is undirected. The cardinality of I i i is called in-degree node v i , and it is indicated with id i while the cardinality of I o i is called out-degree node v i . The Boolean network is dynamic because the configuration c ∈ C is a function of the time. If the status of each node s i (t + 1) is determined only by the value of the states of the nodes of the input neighbor I i i , namely: where each s r represents the state of the node v r ∈ I i i , and this function does not change over time, the Boolean network is called synchronized and the dynamic system is deterministic, in the sense that the successor of a configuration is unique. Equation (1) is given usually as a table of rules T i .
The evolution rule of the Boolean network f : C → C is expressed as the set T of the tables of the local rules T i : namely: Since id i is the cardinality of the surrounding of I i i ,, then the number of rules in the table of rules T i will be just id i and, therefore, the number of rules contained in the rules table of T is: In this way, all the possible combinations of inputs are equal to 2 n t . In the modeling of cognitive processes related to Action-Execution, the connectome with 82 nodes becomes the Boolean network. The nodes of interest are those identified for the chosen behavioral domain, whereby they will have a state of activation as the initial condition. An initial state must be set by assigning to each node a Boolean value. Later, the simulation provides for the transition of the state of each node based on the Boolean rule that is assigned to Boolean variables and the nodes to which it is connected. The state of a node at time t is given by x i (t), and the next state after each iteration is given by www.nature.com/scientificreports/ where x ij are the states of the nodes related to node i . The states of the nodes of the network are updated simultaneously according to this rule. Usually, this process is iterated until the reaching of a stable fixed point or limit cycle. This model allows a systematic exploration of the relationship between brain network structure and dynamics that might otherwise be impossible. In its most simple realization, each node of a Boolean network is connected at random to a set of K input nodes, whereas in our case the connections are provided by the connectome. In addition to the classic methods, further schemes can be introduced based on the idea of a threshold activation. According to this scheme, after choosing an initial state, the Boolean functions associated with the node will not normally be used, we instead calculate the sums of the values of the Input regions. By using this scheme, updating will take place in such a way that if at step t , the ith node is connected to a number of a switched node that is within a given range, then it will be switched on at time t + 1 , otherwise it will be turned off. There are many activation threshold methods, which are suitable for brain and cognitive processes networks modeling. In our case, we will use a consistent organization of the connections, linked to well-known cognitive circuits, as in 59 .

Step 3: local and global cognitive dynamic analysis
The managing of all these variables allowed us to get data on the local interactions among brain areas.
The connections between the various Broadmann areas are given by a connection matrix 82 × 82 (see Table 1) whose values are between 0 and 1. To transform this matrix into an adjacency matrix capable of modeling the system as a not fully connected graph, it is necessary to introduce a threshold x for the values of the matrix, which, below the threshold will be set equal to 0 while above the threshold will be set equal to 1.
Starting from a threshold value of 0.8, it is possible to increase this value by 0.01 every row of the matrix, thus having fluctuating values (which correspond to different values between the nodes).
If we consider the sum of states in a neighborhood σ i = (s i1 (t) + s i2 (t) + · · · + s ik (t)) of the node v i , then we can define the updating function (3) as follows: By varying the threshold of the x matrix and the values of a and b , we can obtain a deep investigation of brain dynamics at the global level of the systems. Once the pattern has been obtained, applying one of the above-mentioned threshold levels, we obtain a connection matrix which highlighted the emerging synchronized circuits. As we will see later, single nodes will exhibit dynamics that can be correlated with the dynamics of other nodes. By calculating these correlations and using a threshold value z , we could identify which dynamics have correlations above that threshold and consequently identify which nodes are correlated with each other. The method allows for the circuits extraction and identification in the network (the nodes which are connected) and in the Brodmann specifications (the involved connected Brodmann areas which have emerged from the dynamics). Many visualizations techniques, at this level, allow us to get the visualization circuit by circuit, and the overall circuit. Analysis of each circuit and of its dynamic behavior were also made. We mentioned that a Boolean network is a graph with N nodes where, at each node, a Boolean variable representing the state of the node (1 or 0; ON or OFF), and a Boolean function that determines the evolution rule associated with it are assigned. To start the simulation, an initial state, assigning to each node a Boolean variable, is chosen. The dynamics provide for the transition of the state of each node based on the Boolean function that is assigned to it, and the Boolean variables of the nodes to which it is connected. So, the system behavior is emergent and unpredictable. As we have said, we will analyze both the local and the global dynamics. In local dynamics, we will investigate the interactions of each node and of its neighbors. In the global dynamics, we will study the dynamics of all the nodes that belong to the behavior under examination. This will allow us to: a. Consider the interactions of each brain region with other regions, discover the rhythms of interaction among the regions, where these exist, or if the interactions are purely random or chaotic. b. Explore the dynamics that lead to the creation of emerging networks that are structured precisely according to the process that is taking place. c. Establish an interaction time, taking care to detect whether emergencies are the result of a latent synchronization process, or whether there are submodules of the brain that synchronize and whether these submodules are in some way coordinated with each other.
Another problem we analyzed in the global dynamics behavior is the stability of the emerged circuits. Obtaining the local dynamic circuits, we try to analyze the global dynamic with an example. In an earlier evolution of BN, we got 10 circuits, considering a = 1 , b = 4 , ( x > 0.9 and a correlation threshold z > 0.87)) as network parameters and threshold values respectively. These ten circuits are considered as a reference configuration to which all the other configurations compare, obtained by changing the parameters a , b and x . In this way, we identified an application in a four-dimensional space R × N × N × R which is the space C of all the circuits. The two thresholds vary in the space R of real numbers, even if the values of the connection weights matrix and the threshold values are conventionally rounded to the second decimal (varying between 0 and 1). Let us denote these discrete spaces regarding the thresholds with H and K . Thus, we have achieved an application f that, takes 4 numbers (integers-a an b , and two numbers between 0 and 1 at the second decimal digit) as input and returns the circuits organization in the four-dimensional space: www.nature.com/scientificreports/ where P = H × N × N × K and C is the set of all possible circuits sets. The problem of global stability is to understand how these circuits vary, at the variation of the 4 parameters. Setting one of the parameters ( a , b , or x equal to a certain value), one of the Ŵ takes an equilibrium configuration Ŵ 0 which we call the cognitive state, assuming that this state is stable. As a , b , or x can be varied, the structure deforms and assumes equilibrium configurations Ŵ 0 ( a , b , x, z ). These are assumed to be discretely dependent on a , b , or x and stable for sufficiently small a , b , x or z . Such a problem can be extensively studied only with a great deal of computational power. In our case, we limited our attention to the understanding of what is going on in the cognitive circuits' stability, taking as example this first set of circuits and comparing all the circuits obtained by varying the parameters to this reference frame, changing one parameter at a time. Our hypothesis is that such circuits can vary slightly within certain limits. So, not all those that vary greatly, can be considered realistic from the point of view of their cognitive reliability. Other features that can be investigated, related to the stability of the circuits, are the concepts of "sensitivity to initial conditions", "damage spreading", and "robustness to perturbations". By changing, mutating or perturbing one node or group of nodes state in the Boolean Network or by mutating their connections, or the rules of the lookup table we can quantify how a casual modification can influence the network, exploring change propagation. After discovering the circuits, we can come back to the real brain and establish whether the emerged circuits correspond to the simulated cognitive functions and whether the shortest path of the network connectome could be considered reliable for detecting how distant brain areas connect together.

Results
Results have been organized along the pipeline of the used methods.
Results of step 1: the structural-functional features of the network. Before starting the simulation for obtaining Boolean Network simulations of the cognitive behavioral model Action-Execution (Fig. 2), the main activation map for this category was accessed at the Brain Map.org database. From the analysis of the metadata for this behavioral category, the Brodmann areas, involved in this process, were extracted and reported in Table 2. These areas, organized in the Sleuth workspace, were visualized with the software Mango, both Brain-Map applications, plotting their results as coordinates on a standard brain in the Talairach space.
The structural features of the network. Standard methods were used to segment a DICOM file and get the connectome. In particular, the following procedure was conducted: Brodmann areas 3, 1, and 2 refers to the primary somatosensory cortex. In the connectome, these areas correspond to nodes 1R and 2R Primary somatosensory cortex (represented by the Brodmann areas 3, 1, and 2 and corresponding to the network nodes 1R and 2R) is involved with the localization of the input stimulus, the evaluation of its intensity, the proprioception and the shape recognition processes. Area 3 receives the information that are then sent to areas 1, 2 and motor areas by the cortico-cortical neurons pathway Brodmann area 4 denotes the primary motor cortex of the human brain. In the connectome, it corresponds to node 8. Located in the rear part of the frontal lobe, the motor cortex is involved in planning, control and execution of voluntary movements of the body, with the function of transmitting to the cells of the nuclei of the cranial nerves and cells of the spinal cord impulses for movements

The parietal cortex The frontal cortex
Brodmann area 5 is part of the parietal cortex of the human brain. In the connectome, it corresponds to nodes 9 and 10. It is involved in somatosensory processing and association Brodmann area 8 is part of the frontal cortex of the human brain. In the connectome, it corresponds to node 15. The area is involved in the management of uncertainty

The anterior portion of the prefrontal cortex The frontal cortex
Brodmann area 10 is the most anterior portion of the prefrontal cortex of the human brain. In the connectome, it corresponds to node 19. Brodmann area 10, the largest area of the human brain and the most unknown, is involved in strategic processes in memory recall and various executive functions Brodmann area 11 is part of the frontal cortex in the human brain. In the connectome, it corresponds to node 21. It is involved in decision making and processing rewards, planning, encoding new information into long-term memory, and reasoning Brodmann area 19 is part of the occipital lobe cortex in the human brain. In the connectome, it corresponds to node 28. In humans with normal sight, extrastriate cortex is a visual association area, with feature-extracting, shape recognition, attentional, and multimodal integrating functions Brodmann area 23 corresponds to some portion of the posterior cingulate cortex. In the connectome, it corresponds to node 35. It communicates with various brain networks simultaneously and is involved in various functions such as human awareness, pain, and episodic memory retrieval

The subgenual area of the cerebral cortex The retrosplenial region of the cerebral cortex
Brodmann area 25 corresponds to the subgenual area of the human cerebral cortex. In the connectome, it corresponds to node 39. Extremely rich in serotonin transporters, this area is involved in vast networks comprising areas like hypothalamus and brain stem, which affect appetite and sleep; the amygdala and insula, which affect the mood and anxiety; the hippocampus, which plays an important role in memory formation; and some parts of the frontal cortex responsible for self-esteem Ectosplenial area 26. It is a cytoarchitecturally defined portion of the retrosplenial region of the cerebral cortex. In the connectome, it corresponds to node 41. This area is bounded externally by the granular retrolimbic area 29 The parietal area of the cerebral cortex www.nature.com/scientificreports/ (1) Acquisition of DTI and high-resolution T1-weighted MRI. Diffusion MRI allows the orientation of fibers in white and gray matter to be determined, exploiting the Brownian motion of water molecules within tissues. This motion can be isotropic-with molecules spreading in similar patterns in all directionsor anisotropic-with molecules following more irregular directional patterns. Isotropic diffusion takes a sphere-like shape, while anisotropic diffusion is ellipsoid-like. (2) Segmentation of the brain's white matter. White Matter Tractography (WMT) includes a set of techniques that allow the connectivity and orientation of brain fibers to be estimated, through the use of images obtained by diffusion MRI. There are numerous algorithms employed, divided into deterministic, probabilistic and global optimization ones. Deterministic algorithms identify a unique trajectory from a selected starting point, and result in a path that connects discrete regions of the brain. (3) Creation of white-matter tractography. Probabilistic algorithms generate considerably more possible trajectories from a single starting point, by connecting them to several voxels and assigning to each trajectory a different probability rate. In the end, global optimization algorithms generate the optimal path between two regions of the brain minimizing a cost function, which generally describes the flow path and its goodness of fit. (4) Segmentation into regions of interest (ROI). The choice of the network nodes is made according to different atlases that allow a segmentation of the brain by different areas to be carried out (Region Of Interest). These atlases encode different areas of the brain related to specific morphological features or functionalities. Different brain segmentations in more or less defined ROIs, will yield networks with 83 or 1000 nodes using the Desikan-Killany Atlas, extended to include the subcortical regions. This is a very subtle choice, according to the aim of the analysis, as indicated in 4 . Xia et al. 72 . The lack of a universally accepted pattern for the segmentation is a well-known problem. For the purposes of this study, we chose to use 82-node networks, reported in Xia et al. 72 . (5) Construction adjacency matrix to determine the connection between the various structures. Anisotropic diffusion analysis is carried out using various methods. The most used are DTI (Diffusion Tensor Imaging) and DSI (Diffusion Spectrum Imaging). Measuring diffusivity involves the introduction of a gradient, or a variation of the magnetic field in different directions. The motion of molecules then causes changes in the intensity of the emitted signal. The gradients are applied both before and after the pulse at 180°. DTI spreading at a given point (voxel) within the brain is quantitatively described by a diffusion tensor. At least six measurements, calculated from a 3 × 3 symmetric matrix with six degrees of freedom, are required to construct this tensor. A connectome mapper allows the 5 steps of the procedure. It uses FREESURFER for the segmentation, FSL to display dicom files, DTK TRACKVIS to determine the tensor, and to draw the fibers and CONNECTOME to create the matrix. In the fragmentation phase, it would be possible to use the following algorithms to calculate the adjacency matrix: Brodmann area 39 is part of the parietal cortex in the human brain.
In the connectome, it corresponds to nodes 63 and 64. Damage to Brodmann area 39 plays a role in semantic aphasia Brodmann area 45 is part of the frontal cortex in the human brain. In the connectome, it corresponds to node 75. Together with Brodmann area 44, it comprises Broca's area, a region that is active in semantic tasks, such as semantic decision tasks (determining whether a word represents an abstract or a concrete entity) and generation tasks (generating a verb associated with a noun)

The orbital region of the frontal cortex
Brodmann area 47 is part of the frontal cortex in the human brain.
In the connectome, it corresponds to node 79. Brodmann area 47 is involved in the processing of syntax in oral and sign languages, musical syntax, and semantic aspects of language From a study previously carried out 74 , one sample DICOM file, from a healthy control patient, was used. The segmentations were conducted on T1 sequences, acquired by a 3 T GE, with a voxel resolution of 1 × 1 × 0.5 mm. Each slice had a 256 × 256 resolution, while total slices were 368. For the tractography process and the fibers extraction, DTI sequences were used, with a gradient number of 27 (the higher the number, the better is the resolution). The resolution of the voxel was 1 × 1 × 1. The resolution of the slice was again 256 × 256. There were 80 slices. The brain connectome was visualized by a multiplatform MATLAB (The MathWorks Inc., Natick, MA, US) toolbox, called BrainNet Viewer 72 , with an easy-to-use Graphical User Interface (GUI), which allows for the use of a visualization platform and the generation of a brain connectome with a segmentation Atlas, organized according to the Brodmann areas. Deriving many functions from the Statistical Parametric Mapping 8 application (SPM, http:// www. fil. ion. ucl. ac. uk/ spm/), this toolbox allows for the downloading of brain surface, node, edge and volume files. By directly downloading the '.pial' file of hemisphere mesh, generated using FreeSurfer and the '.mesh' files, generated using BrainVISA (http:// brain visa. info/) 75 , a connectome can be built. This tool has been chosen for the reason that it has been modeled on the automate Talairach registration procedure, developed by the Montreal Neurological Institute 76,77 , which allows the correlation between the individual volume and an average volume composed of a huge number of aligned brains to be maximized.
The functional features of the network. The transformation matrix into the Talairach space was used to define exactly the ROI related to the Action-Execution points of activation in the functional maps. After these processes, the brain connectome and its connection matrix were processed in the Mathematica environment (Fig. 3) 78 . Table 2 shows the Brodmann areas corresponding the ICAs behavioral domain area of the connectome visualized within the Talairach atlas 79 by the Mango application software. Figure 4 shows the Action-Execution behavior from the Brodmann atlas and the corresponding node networks in the connectome.

Results of step 2: simulation runs.
In starting the simulation, we decided to activate only the 16 Brodmann areas belonging to the chosen cognitive process. Therefore, the initial condition of the Boolean network was as follows: Parameters of the network are as follows: 1 ≤ a ≤ 13 , 1 ≤ b ≤ 13 , steps of the simulation = 500.

Results of step 3: the emergence of brain circuits.
To highlight the complex emerging dynamics and the presence of structures in the simulations carried out, a model simulation ( b = 4 , network threshold x > 0.9 and, correlation threshold z > 0.87, same initial conditions), in order to consider a sufficiently complex but not chaotic behavior, is analyzed in depth. In this simulation, nodes always turned on and off have been removed, as they do not take part in the dynamics that is going on if not in a subtle way, which we will analyze in the next paragraphs. After the irrelevant nodes have been deleted, we sought correlations among those that changed over time, by means of a correlation matrix (for b = 4), as shown in Fig. 7.
The correlation matrix provides us with a measure of how the performance of a brain area is related to the performance of another area. This value is between 0 and 1. When the correlation between the two areas is equal to 1, the two areas are perfectly synchronized. If we insert a z threshold, selecting the pairs of areas that have a correlation higher than z , then we can identify groups of areas that have fairly similar behavior and we call these areas brain circuits. The model simulation gave 10 emerging circuits. The number of areas for each of these circuits is as follows: {2, 2, 2, 2, 2, 2, 2, 2, 3, 3}. Each circuit connects with up to three nodes of the brain connectome. In Table 4, the emerged circuits are reported, while in Fig. 8, the representation of the emerged circuits on the brain connectome is shown. Table 3. Temporal emergence of periodic behavior in the brain connectome.    Table 5, results about the reliability of the emerged circuits agreeing to the cognitive behavior Action-Execution is evaluated, considering the ROIs normally used in the previously-mentioned real process.

Circuit analysis
Let us investigate whether there is any relationship among the emerged circuits for this simulation. The easiest and quickest way is to adopt the intersections between the areas of the emerged circuits. In Table 6 we show the intersection between these 10 circuits.
A different way to represent this information is in the form of an intersecting array as follows (Fig. 9). Given two circuits r and s , we can calculate the ratio between the number of elements i rs contained in the intersection between r and s , divided by the number n r of elements of r . For example, i 2,7 = 1 2 , while i 2,9 = 2 2 = 1. If we carefully analyze these circuits, we notice some interesting processes:   The emergence of circuits. Using the simulation chosen as an example, and starting from the 16 active brain areas identified for the Action-Execution cognitive process, taken as initial data, let us analyze how the circuits emerge along the evolution of the spatio-temporal structure. After 250 time steps, the system evolves with a periodic configuration (Fig. 10) Of these, only 2 {24L, 35L} and {R 42, R 45} are in common with the final circuits. We may think that these circuits, once they are formed remain even later. In reality it is not so. In fact, if we analyze the spatial and temporal pattern between 50 and 100, we observe the formation of more than 16   www.nature.com/scientificreports/ previous example. Also with the parameters a = 1 and b = 4, we only consider the 16 active areas initially. After just 16 steps the system reaches the periodic attractor. The oscillation period is 6 time steps. Let us now analyze the circuit of Action-Execution as a distinctive one, as represented in Fig. 11, allowing its evolution at the local level of the BN, the realization of the connection matrix, the detection of its main circuits and lastly, their representation on the brain connectome.
To obtain the emerging circuits of this simulation, it was necessary to slow down the threshold of the connectivity ( z = 0.6 ), as a very high level of connections may hide the process (see Table 7).
The dynamics of the emerging circuits: stationary, oscillatory and rhythmical behavior. To better understand the nature of these emergent circuits, let us examine another case in which a = 1 e b = 6 maintaining the same threshold x = 0.9 . Results about the evolution of the Action-Execution circuit and the relationships that each area engages with the other ones by means of the phenomenon of binding together by symmetry 80 are presented in what follows. We have identified oscillatory, stationary and rhythmical repetition of dynamical behavior for the emerging networks. Nodes always switched on or off are as follows : {1, 3, 5, 9, 12, 15, 16, 21, 24, 35, 54, 56, 65, 73, 74, 78, 82}.
The connection matrix for those nodes that vary in time is represented in Fig. 12, together with the connection matrix for circuits.
The emergence of circuits at the global level. By modeling the brain dynamics with Boolean networks, we detected the emergence of cognitive circuits, which are very large-scale-networks. From the simulations carried out extensively for all system parameters, we found many circuits. Usually, the circuits emerge quickly, in large numbers. For the chosen reference, we found that 9 circuits emerged from 0 to 50 simulation steps. Of these, two are in common with the circuits detected at the end of the simulation. 16 circuits emerged between 50 and 100 simulation steps, but only 1 is in common with the circuits detected at the end of the simulation. Subsets of circuits that emerged at the end of the simulation are related to the huge circuits emerged at the beginning of the simulation.
Let us try to analyze the results on the BN global dynamics behavior using the same simulation as an example. Considering a = 1 , b = 4 , ( x > 0.9) as network parameters and threshold values respectively, we got 10 circuits (Fig. 13), considered as a reference configuration to which to compare all the other configurations, obtained by changing the parameters a , b and x . This circuit configuration was chosen because it is very close to the physical organization of the behavioral map, drawn from the ICAs standard. It could be represented also in the brain and it is realistic from a cognitive point of view as a computational model.
In our case, we limited our attention to the understanding of what is going on in the cognitive circuits stability, taking this first circuit configuration as an example and comparing all the circuits obtained by varying the parameters to this reference frame, changing one parameter at a time.
Our hypothesis is that such circuits can vary slightly within certain limits. So, not all those that vary greatly, can be considered believable from the point of view of their cognitive reliability.

Influence of the connection threshold on the stability of the emerging circuits
Let us change the Boolean Network threshold for the following values of the x parameter: Unlike the case considered as reference with x = 0.90 , for which we obtained 10 circuits, for x = 0.89, 15 circuits emerged from the simulation run. To evaluate the possible overlapping among the two different simulation runs, and consequently among the two sets of emerging circuits, an intersection between all the possible circuits of the two simulations has been done. We can represent these intersections using the same method of Fig. 9, obtaining a matrix shown in Fig. 14.
As can be seen, there is no relationship between the two groups of circuits. However, it is interesting to observe that, as occurred at the local level, a new circuit, composed of two areas {24L, 43L}, appears. These two areas already had an important role in the simulation, taken as a reference model. In fact, the new emerging circuit has at least one element in common with six circuits of the reference case. Instead, considering a threshold level of x = 0.91 , we observe the emergence of 22 circuits, of length: {2, 2, 2, 2, 2, 3, 4, 5, 6, 6, 7, 9, 11, 12, 12, 13, 13, 15, 15, 16, 17, 19}. Therefore, in addition to the presence of small circuits made up of at least of three elements, we discovered the presence of larger circuits, composed of many elements. The relationships between these new set of circuits and the circuits emerged for the reference model simulation, are graphically expressed as follows (Fig. 15).
Circuit 1 (repeated 3 times) and circuit 3 (contained in 5 of the new circuits), continue to exist as sub-circuits of larger circuits. So, the threshold according to which the Boolean network is built changes the topology of the   Table 7. Circuits emerging from the Action-Execution evolution with a BN. We hypothesize that these circuits represent the firing over time of all areas involved in the simulated process, with oscillations between areas.

Connectome nodes
Brodmann areas Correlation threshold. The correlation matrices are used to identify circuits, analyzing the related areas.
In particular, an area is related to another area if the correlation of the signals associated with the two areas is greater than a certain threshold z . Given an area r , and a threshold z , the set of all areas A r whose correlation with the r area is greater than z represents a circuit, adding r . Each circuit is counted only once. It is important to note that two areas a 1 and a 2 can belong to the same circuit, even if the correlation between them is less than z . This is the case of the reference circuits 2, 7 and 9. Circuits 2 and 7 have two areas correlated with z > 0.87 . Their union represents a new circuit (circuit 9) that contains three elements. Two pairs, corresponding to circuit 2 and 7 have a higher correlation of z , while the third pair has a correlation < z . As z decreases, the starting circuits are maintained, and maybe are incorporated in larger circuits, while for growing values z , these circuits instead tend to clear. The new circuits will be subsets of the starting circuits. Let us analyze in detail different correlation thresholds. For z = 0.86 , there are still 10 circuits and obviously identical to those of the reference case. The case where z = 0.85 is more interesting where there are now 15 circuits of length {2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 6}. It was possible to observe that the first circuit continued to subsist; circuit 2 continued to exist individually and was also part of a larger circuit. Even circuit 3 became part of a larger circuit. Circuit 4 continued to maintain its independence, while 5 appeared as a single entity and also part of a larger circuit. Let us observe, however, how the circuits are developed in the direction of increasing thresholds. If z = 0.88 , there are now 8 circuits. All the reference circuits continue to survive, with the exception of circuit 5, which disappears, and circuit 10 that has lost one of the three areas coinciding with circuit 6. If we further increase the threshold to 0.89, we observe that now there are now 5 circuits as seen from the Figure xxu1. The only circuits that tend to survive are circuits 3, 4, 6, 7 and 8. While circuits 9 and 10 lost each area of the three they had previously. The result did not change until value z = 0.92; while for z = 0.93 the only circuit that continued to survive was 8, corresponding to the two Brodmann areas {42R, 45R}.

Discussion and conclusions
To model structural and functional aspects of brain dynamics better and to be able to grasp cognitive behavior, we used Boolean Networks for the following reasons: (a) The attractiveness and ease-of-use of the BN model 81 which, by simplifying the complexity of brain dynamics, allows the emergence of circuits and chain of circuits (the network's attractors ) which in turn represent stable (or metastable) sets of states of the brain toward which the system converges. Therefore, evolving toward attractors means to reach a specific behavioral goal (probably by a completely mechanical internal dynamics). This approach helps to predict the brain's organization at different scales (the single circuit dynamics at the level of each brain area, the emergence of circuits that carry out the specific functional task, the completely new emergence of circuits, which could shed light on the long-range informational process in the brain). In this way, we have a running idea of what is a cognitive-functional attractor and what is the structural side of the system that allows its dynamics.   3,82,83 . Therefore, comparing simulated dynamics with observed brain networks of a certain cognitive behavior may help in evaluating their reliability in predicting and showing missed information which could not be seen, by using real life methods. BN helps in predicting the process, picking up the qualitative and quantitative changes in the dynamic expression of the cognitive behavior. (c) This model allows observation of the structure and lengths of the behavior, especially when expressed in recurrent patterns of activation. This in turn allows the determination of the basins of attraction of the brain circuitry, and to analyze the flow of information through the brain, thus determining the phase transition between chaotic, transient, ordered patterns, also discovering the edge of chaos [84][85][86] , which is where complex computation occurs. In addition, topology is of primary importance, owing to the fact that different dynamic behavior and flow of information can occur according to different parameter values, initial conditions, thus giving account of the brain as a dynamic non-linear 87 entity. (d) BN well fits with the anatomy of the brain, represented as a connectome, for understanding basic brain functioning, and, therefore, both the structural and functional sides of the brain can be properly replicated in a biological-like structure. The finite set of the activation rules and updating scheme of the BN, together with the possibility of varying the synchronous or asynchronous method, and the initial conditions allow a reduction of the potentially infinite number of processes a brain can carry out. For each node, the activation rules are determined by the activation states of the other nodes in the previous transition step, synchronously or asynchronously. These critical features of the BN dramatically influence the number and characteristics of the attractors. As brain areas activation rules and neighboring interaction among areas are not known, by using BN different activation rules and schemes can be tested, which best fit with conditions under which the simulated network would be able to capture the observed cognitive processes.
Another point of discussion is the choice of different parameter values of the NB, especially for threshold values. The threshold became the filter to see how nodes of the brain are individually correlated with the other nodes. It was crucial to establish the parameters according to careful neurocognitive modeling. This process would allow understanding the relationships between brain areas that are activated for the cognitive process we are describing. This, however, would require that it can be known, a priori, what these connections are from the neurophysiological point of view. As in most cases there is no data in this context, we tried to do some evaluations of the best possible parameters. We can state that the rule of evolution of the Boolean network is crucial to the stability of the brain circuits, which, therefore, significantly change their characteristics. For increasing values of b and decreasing values of b , we observed the declining of many circuits, while for growing values of b many circuits disappeared, although a few of them remained as subcircuits of larger circuits. For the simulations we realized, putting the threshold at 0.9, the system produced no important interactions, it did not reveal any This is a very simple circuit with a large number of nodes which rythmically turn on and off at each step of the simulation Table 8. Emerging circuits and their dynamical behavior with a = 1 e b = 6 maintaining the same threshold x = 0.9.
Scientific Reports | (2022) 12:16543 | https://doi.org/10.1038/s41598-022-20979-x www.nature.com/scientificreports/ meaningful structure. This meant that, the nodes of the brain circuitry having very strong connections, did not imply having effective interactions or correlations that can emerge. Placing the correlation threshold at 0.80, interaction among areas and many brain circuits emerged. We chose to put the correlation threshold at 0.87 in order to get a rich variety of connections. This question is strictly related to the brain local and global dynamics, and the observation of phenomena at different scales. According to this issue, we noticed mixed phenomena. The effects of the threshold levels changed the results of the simulations. In fact, if a very low threshold level allowed the emergence of the circuits at the local level, at the global level, this has not happened. This differs from the  www.nature.com/scientificreports/ previously-argued hypothesis at the global level, for which we have shown that the average low threshold reveals an acceptable number of circuits. On the contrary, the high threshold increases the number of circuits. The phenomenon is certainly linked to the scale of observation of the dynamics that are emerging. Therefore, it can be stated that, according to the level of scale, the functions of the thresholds can vary, providing results of completely different brain circuital connections (although they can be related). So, at the global level, a low threshold can make the circuits disappear, while at the local one, they continue to function. From the mathematical point of view, we are dealing with two different measures, in respect of which, there is no correlation. If we were to use the global threshold for the local level, there would be no local circuits. Lowering the local threshold is as if we accept the existence of local weak signals and strong global signals. It is a modeling outcome. We do not know if it also has a corresponding phenomenon in neuronal areas at different scales. In fact, the correlation that we use is a coarse modeling, interesting at the global level, but less significant at the local one.
In conclusion, we can say that modeling with Boolean Networks makes sense from the neurocognitive perspective. In fact, in the cognitive process we have simulated, it is known that the first circuit affects the somatosensory cortex and the entorhinal cortex. This means that the two areas are activated for bringing spatial, visual and tactile information to the cognitive system (from the somatosensory cortex), which relates in some way to the stimulus localization. Then stimuli are conveyed to the entorhinal cortex system, where the grid cell systems provide information on the agent's localization in the space in which he/she is in, especially thanks to the head position sensors, in correlation with a series of visual and auditory stimuli. These data, in turn, allow information to be obtained on the shape and size of objects on which he/she has to actually carry out the action. Furthermore, the entorhinal cortex connects with other systems, including the hippocampus, which plays a significant role for storing processes and actions that have already been realized. Brodmann area 43 L is common to circuits 2, 7 and 9, while Brodmann area 24L is common to circuits 5, 6 and 10. This might suggest that there are 6 circuit clusters of a higher level: the first four coincide with circuits 1, 3, 4 and 8, connected to Brodmann areas which do not communicate with other circuits. Two cerebellar areas are made first, from circuits 2, 7, 9, while the second of 5, 6, 10 circuits. In reality, circuit 9 has the two sub-circuits 2 and 7 whose areas are strongly correlated with each other via area 43L, while the other 2 areas are weakly correlated. Similarly, for circuit 10, which contains both circuits 5 and 6, weakly connected between each other and strongly by Brodmann area 24L.
A complex landscape of how cognitive processes emerge is given by the phenomenon we interpreted as the emergence of metastable structures, which gave computational data on large-scale networks emergence and on www.nature.com/scientificreports/ how the brain creates structures of remote communication. This is a metastable structure 88 , with a hierarchical organization, where each level allows for the emergence of brain organizations which behave at the next superior level.
In the dynamics of evolution for this simulation, we observed both stationary and oscillatory behavior, and repetitions of rhythms that can be ascribed to the different ways brain areas have to interact, in the Action-Execution simulation process. At the global level, after evaluating the whole period of the network, we denote that the above-described dynamics achieve their behavior each 48 steps of simulation, with complex organization at the edge of chaos. This problem can be extensively studied only with a great deal of computational power. In our case, we limited our attention to the understanding of what's going on in the cognitive circuit stability, taking as example this first circuit configuration and comparing all the circuits obtained by varying the parameters to this reference frame, changing one parameter at a time.

Data availability
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.