An Evolutionary Game Theory Model of Spontaneous Brain Functioning

Our brain is a complex system of interconnected regions spontaneously organized into distinct networks. The integration of information between and within these networks is a continuous process that can be observed even when the brain is at rest, i.e. not engaged in any particular task. Moreover, such spontaneous dynamics show predictive value over individual cognitive profile and constitute a potential marker in neurological and psychiatric conditions, making its understanding of fundamental importance in modern neuroscience. Here we present a theoretical and mathematical model based on an extension of evolutionary game theory on networks (EGN), able to capture brain's interregional dynamics by balancing emulative and non-emulative attitudes among brain regions. This results in the net behavior of nodes composing resting-state networks identified using functional magnetic resonance imaging (fMRI), determining their moment-to-moment level of activation and inhibition as expressed by positive and negative shifts in BOLD fMRI signal. By spontaneously generating low-frequency oscillatory behaviors, the EGN model is able to mimic functional connectivity dynamics, approximate fMRI time series on the basis of initial subset of available data, as well as simulate the impact of network lesions and provide evidence of compensation mechanisms across networks. Results suggest evolutionary game theory on networks as a new potential framework for the understanding of human brain network dynamics.

MRI data was acquired on a 1.5 T (Philips Intera) whole-body scanner. Resting-state fMRI data included 178 volumes with 33 axial slices covering the whole brain, acquired via a T2 BOLD-sensitive multi-slice echo planar imaging (EPI) sequence (T R/T E = 2.5s/32 ms; field of view= 22 cm; image matrix = 64 × 64; voxel size = 3.44.44 × 3.8 mm 3 ; flip angle = 75 • ). Structural imaging was performed using a whole brain T1-weighted Fast Field Echo 1 mm 3 sequence (T R/T E = 30/4.6 ms, field of view = 250 mm, matrix 256 × 256, flip angle = 30 • , slice number = 150). T2-weighted Fluid Attenuated Inverse Recovery Images (FLAIR) were also acquired to assess participants white matter integrity. Participants were provided with earplugs. Particular care was taken to minimize head motion via vacuum cushions and custom made padding. Both structural and fMRI data have been carefully checked for data quality.

fMRI data preprocessing
Functional image preprocessing and statistical analyses were carried out using SPM8 software (Statistical Parametric Mapping; http://www.fil.ion.ucl.ac.uk/spm/) and MATLAB 7.5 (Math-Works, MA, USA). The first five volumes of functional images were discarded for each subject to allow for steady-state magnetization. EPI images were slice-time corrected using the interleaved descending acquisition criteria, and realigned and re-sliced to correct for head motion using a mean functional volume derived from the overall fMRI scans. Subject whose head motion exceeded 1.0mm or rotation exceeded 1.0 • during scanning were excluded. In order to obtain the better estimation of brain tissues maps, we implemented an optimized segmentation and normalization process using DARTEL (Diffeomorphic Anatomical Registration using Exponentiated Lie Algebra) [1] module for SPM8. Briefly, this approach is based on the creation of a customized anatomical template built directly from participants T1-weighted images instead of the canonical one provided with SPM (MNI template, ICBM 152, Montreal Neurological Institute). This allows a finer normalization into standard space and consequently avoids under -or overestimation of brain regions volume possibly induced by the adoption of an external template. Hidden Markov Random Field model was applied in all segmentation processes in order to remove isolated voxels [2]. Customized tissue prior images and T1-weighted template were smoothed using an 8mm full-width at half-maximum (FWHM) isotropic Gaussian kernel. Functional images were consequently non-linearly normalized to standard space and a voxel resampling to (isotropic) 3 × 3 × 3 mm were applied. Linear trends were removed to reduce the influence of the rising temperature of the MRI scanner and all functional volumes were band pass filtered at (0.01Hz < f < 0.08Hz) to reduce low-frequency drift. Finally, a CompCor algorithm has been applied in order to control physiological high-frequency respiratory and cardiac noise [3].

EGN-B Model
The model proposed in this paper assumes that each agent is an aggregated (brain area) of interacting elementary components (neurons). Each area is assumed to be able to make decisions, and it is characterized by its activity level or state. The state, belonging to the set [0, 1], quantifies the activity level of the area at a particular time. Dynamically, each area compares its activity level with others and it changes its state accordingly, in order to maximize a certain payoff function based on activation and inactivation strategy mechanisms defined below. The above mechanism is naturally embedded into the replicator equation on graphs [4], describing the EGN-B model as follows: where: p A v,w and p I v,w represent the payoff of activation and inactivation, respectively. Moreover, α v,w and ι v,w are the activation and inactivation propensities of payer v with respect to player w. Finally, x v ∈ [0, 1] is the level of activation of the brain area v.

Estimation method
In order to characterize the system (1), we must estimate 3N 2 parameters, namely a v,w , α v,w and ι v,w , using the available observed data (the fMRI recordings correspond to the behavior of particular brain areas). The estimation procedure is constrained: in fact, the entries of adjacency matrix a v,w must be non-negative, while the propensities α v,w and ι v,w must have the same sign for each ordered couple (v, w). In this work we relax these constraints in order to tackle with an unconstrained optimization problem. To this aim, we set α v,w = ι v,w ; in this way, the threshold d v is set to 0.5. Moreover, since the parameters a v,w quantify the strength of connections, then we assume that α v,w and ι v,w belong to the set {−1, +1}. Moreover, parameters α v,w and ι v,w in (1) are always multiplied by the non-negative parameters a v,w , thus affecting the identifiability of the model. For these reasons, we can instead estimate the following N 2 parameters: Since α v,w = ι v,w ∈ {−1, 1}, the parameters reads as: The result of the estimation method is represented by the matrix A = {a v,w } ∈ R N ×N which simultaneously contains two important informations: • since A = |A |, it encapsulates the network structure; Notice that, if a v,w = 0, then w does not influence v. On the contrary, if a v,w = 0, then |a v,w | denotes the strength of their interaction, while sign(a v,w ) denotes if both activation and inhibition propensities are positive or negative.
Starting from these premises, we can rewrite the EGN-B in (1) as follows: Last equation shows that the system (2) linearly depends on the parameters of the signed EGN-B connectivity matrix A , and it can be rewritten in an alternative form: where Moreover, the model (3) can be discretized using the Euler method with time step τ S : ) is the approximation of the time derivativeẋ v , and τ S is the sampling time.
According to (4), the parameters of the model (i.e. the entries of the signed EGN-B connectivity matrix A ) can be estimated using a least square optimization approach. In particular, we indicate with z v (t k ) the observed fMRI data collected with sampling time τ S for T time instants. By substituting x v (t k ) with z v (t k ) for each k ∈ {0, 1, . . . , T − 1}, we can introduce the vector: . . .
thus obtaining the following matricial form: where θ = [θ 1 , θ 2 , . . . , θ N ] . Finally, we can estimate the parameters θ by solving the following problem:θ Functional in equation (6) is convex and, under certain hypotheses (see Technical issue and remarks paragraph at the end of this section), it has only one solution: whereθ are the estimated parameters based on the observation z.
Notice that the state variables of the dynamical model involved in this work are bounded in the set [0, 1]. For this reason, data before interpolation have been normalized in the set [0.3, 0.7].
Technical issue and remarks Notice that Y (z) ∈ R N T and U(z) ∈ R N T ×N 2 . Then (U(z) U(z)) −1 is well-defined if rank(U(z)) = N 2 . For this reason, it is a good practice to set T > N , thus allowing U(z) to have more rows than columns; this increases the probability to deal with a well-posed problem. Unfortunately, the number of available measurements is often small; fMRI data scan makes sense if the subject remains in the identical position for the whole duration of the experiment, that can persist in the interval [5,10] minutes at a sampling time in the range [2,3] seconds. Anyway, the bandwidth of fMRI signals is B = 0.1 Hz and the measurements in this work are sampled with a frequency 1 Hz. Thanks to the Nyquist-Shannon theorem, since the f S ≥ 2B all the frequency information of signals are kept after the sampling process. Hence we can interpolate the available data in order to have a finer time discretization, thus increasing the value of T up to the desired target. In this work, we have interpolated the fMRI signals by using splines, adding 100 samples between each couple of time consecutive real data.
We remark that the dimensionality of the problem can become a problem as long as the size of the considered network increases [5]. In order to avoid identification problems when dealing with huge networks, one can imagine a multiple steps identification where the weak connections are removed in order to reduce the size of the problem using a priori information like the knowledge of physical and anatomical connectivity.
However, a sufficiently large number of nodes is required to guarantee satisfactory model performances. To show this, we considered several different size models, increasing from 2 to 77 nodes. This investigation has been performed for all the subjects involved in the study. Thereafter, for each estimated model, we computed the average fitting error as a function of the size of the network, as depicted in Figure S1-A. Notably, as long as network size increases, the fitting error decreases, thus improving the model performances. Moreover, the cross correlation between real and simulated data for each network size has been evaluated. In Figure S1-B, we report the minimum (cyan) and the maximum (magenta) values of the cross-correlation curves. We observe that the min/max cross-correlation range increases as the number of the network nodes increases. Furthermore, when the network size is sufficiently large (> 65 nodes), both minimum and maximum cross-correlation values are positive, thus combining the high fitting model performances with a strong similarity between simulated and real data. Finally, for more than 70 nodes, the cross-correlation raises to 1, making reasonable the extraction of the 77 nodes network from the atlas performed in this study.

Comparison with linear systems
In order to evaluate the performances of the proposed EGN-B model, we performed a comparison with standard linear systems. Firstly, for each subject in the available dataset we estimated the EGN-B adjacency matrix using a number of estimation samples N E lower than the number of samples of the full recording of a single subject (240 samples). Thereafter, we simulated the model using the estimated adjacency matrix for 240 time instants, and we calculated the mean square error between real and simulated data. This error quantifies how much our model is Figure S1: Performance of the model for different network sizes. For all subjects involved in the study, the adjacency matrix of the model has been estimated by using an increasing number of nodes (from 2 to 77). In (A), the average fitting error between simulated and real data is shown, while in (B), the minimum (cyan) and maximum (magenta) cross-correlation curves between real and simulated data are reported.
able to predict the brain dynamics using a reduced number of samples N E . For each size of the estimation set N E , we evaluated the mean, standard deviation, best and worst cases obtained within the group of available subjects under investigation. Moreover, we counted the number of outliers N O within this groups, as the number of subjects whose prediction error is bigger than the mean plus 3 times the standard deviation.
Afterwards, we performed a similar procedure based on a linear system, defined as: This time, we estimated the dynamical matrix M using the same approach employed for the signed EGN-B connectivity matrix A .
The results obtained using the two models are reported in Table S1. Dark blue lines refer to the results obtained with the EGN-B model, while light blue lines reports the results for the linear model. The mean error is always lower when the EGN-B is used to predict the brain dynamics of the subjects (except for the case with N E = 150). The mean obtained using the EGN-B model does not present strong differences when N E varies -its order of magnitude is always 10 −2 -as it happens for the linear model. Moreover, the results obtained using the proposed model are more robust than the linear case: indeed the latter shows highly variable performances compared to the the EGN-B model (see the standard deviation and worst case columns of Table S1). Finally, the number of outliers in the investigated groups is always null for the EGN-B case when the estimation set is small (N E ≤ 140) compared with the linear case when there is always at least one outlier.  Table S1: Comparison between EGN-B and linear model based on prediction error. For each size of the estimation set N E , we report the mean, standard deviation (Std), the worst case (WC), the best case (BC) and the number of outliers N O within the group of available subjects under investigation. The results obtained using the EGN-B and the linear model are reported in the dark blue lines and light blue lines, respectively. Figure S2: Functional connectivity (correlation) simulation and temporal prediction. Functional connectivity matrices built using an incremental amount of available data points are shown (1/3, 2/3, 3/3), leading to a matrix almost identical to the original one when 100% of the data are modeled (3/3). Interestingly, the application of EGN-B on the entire set of time points also allows to simulate data beyond the available dataset, as in the case of the connectivity matrix on the far right (4/3), showing changes in future network behavior by simulating an additional 33% of data points.