Functional brain networks are networks created with a help of fMRI measurements of the in vivo brain activity [3], [4]. I elaborated a dataset, which contains functional brain networks of young participants, healthy elderly participants and elderly participants with diagnosed Alzheimer disease. All networks were measured at the three different correlation thresholds. In this paper I present a data driven mathematical model of functional brain networks. It is based on the threshold related shape of the degree distribution. The model is numerically simulated and results of the simulation are compared to the real dataset.
Functional Magnetic Resonance Imaging (fMRI) is a technique for gaining high resolution images of neural activity in the brain [3], [4]). FMRI images are captured in a series of two dimensional slices, with each slice representing a cross section of the brain less than 10 mm thick. A single slice is comprised of a rectangular grid of discrete 3D regions (3 × 3 × 10 mm) known as voxels (volumetric pixels). A full 3D image of the brain is achieved by combining these slices together. fMRI is an ideal technique for deriving functional connectivity. One can ask to which extend spatially distinct regions of the brain exhibit similar behavior over time. By modeling this functional connectivity as a network, we can explore the ways in which regions of the brain interact, and use techniques from the graph theory to evaluate the topological characteristics of these functional networks of interaction.
In this paper I used the brain fMRI data collected by Buckner [1], (data set no. 2-2000-118W from the fMRI Data Center: http://www.fmridc.org). The participants were divided into three groups: healthy young (HY) participants, healthy elderly (HE) participants and elderly participants with diagnosed Alzheimer disease (AD) (AE group). Structural and functional MRI data were acquired from 41 subjects in total. The HY group had 14 subjects (9 females/5 males) with the mean age 21.1 years (SD 2.0). The HE had 15 subjects (9 females/6 males) of the mean age 75.1 years (SD 6.9). The AE group had 12 subjects (7 females/5 males) of the mean age 77.1 years (SD 5.3). There was no statistically significant difference in the mean age of the latter two groups.
The standardized data were then used in [3], [4] to create functional brain networks for each participant in all three groups. These data were further elaborated by me.
Functional brain networks, contrary to the structural
neuronal brain networks, are temporal networks. Certain
type of the functional brain network exists only during
that time, when the brain is involved in the cognitive task
and reflects the functional cooperation of different brain
areas. Since the smallest unit of the measured fMRI
signal is an integrated signal of the neurons contained in one
voxel, voxels are thus natural candidates for the nodes of
the functional brain network. If the two voxels
functionally cooperate (based on the underlying physical
connectivity), the measured signal is highly correlated over time.
To measure the amount of the signal correlation, the
Pearson correlation coefficient is calculated for all the voxel
pairs:
r(i, j) =
(
Simple measures, that characterize the network in general are averages: such as average degree, density, average shortest path, average clustering coefficient, etc. Usually, these simple measures are not sufficient and one have to rely on distributions, such as degree distribution for example, to acquire more detailed network properties. 2
In this section I complete a basic network analysis of the measured data, already done and published by [3], [4]. The networks are constructed with respect to the three different correlation thresholds, namely θ1 = 0.819398, θ2 = 0.899876 and θ3 = 0.962249. McCarthy et al. [3], [4] have calculated average properties, i.e. number of Degree distribution HY, threshold 0.899876 gamma=-1.1407 nodes, density, degree, clustering coefficient, path length, local and global efficiency, small world index, assortativity for all classes of functional networks. These values were compared across each of the three groups, in order to find differences related to the age and the presence of AD. I concentrated my attention to the degree distributions in all three groups of participants at each of the three thresholds. The reason is, that the dynamical functional brain network model is based on this.
First, the whole brain networks of the healthy young
(HY) participants are described. Then I mention also the
other participant groups. At the beginning I have to state,
that network is scale free if it has a power law degree
distribution (
For the lowest correlation threshold θ1 = 0.819398 one can see in Fig. 1, that the degree distribution does not have a power law character. Thus, the functional brain networks are not scale free and the tail of the distributions is not long enough to estimate the power law scaling exponents correctly.
Degree distribution HY, threshold 0.818398 0 -2
For the correlation threshold θ2 = 0.899876 the degree distribution reveals more pronounced power law tail , with 2 the average scaling exponent γHY = −1.14 (see Fig.2). The scaling exponent is calculated as an average of all scaling exponents of all distributions for the threshold in question.
The situation changes dramatically for the highest correlation threshold θ3 = 0.962249. The functional brain networks now have a well defined scale free structure, reflected in the power law degree distribution with the average scaling exponent γ = −1.36 (see Fig.3) with the individual differences in the interval [−1.7812, −0.9346].
The same statistical distributions as for the previous group were analyzed for the HE group of participants. For the lowest correlation threshold θ1 the degree distribution shows similar features as in the group of the healthy young participants, with the exception, that the individual differences are more pronounced.
For the medium threshold θ2 the degree distribution reveals more pronounced tail in the log–log plot and shows more variability in individuals then the similar degree distribution of the HY group. Statistical analysis of the networks generated for the highest threshold θ3 shows that the degree distribution seems to be scale free, but with more individual differences than in the HY group. Average scaling exponent is γHE = −1.3609 and all individual scaling exponents are in the interval [−2.0396, −1.0500].
We have also analyzed the functional brain networks of the elderly people with diagnosed mild or very mild Alzheimer disease (AE group) for all of the three thresholds. In comparison to the first two groups, namely HY and HE, we have noticed greater individual differences. Even for the highest threshold not all of the networks are scale free.
For the θ1 correlation threshold, the networks are not scale free, There are more pronounced power law tails of the degree distributions at the θ2 correlation threshold and the average scaling exponent is γ2 = −0.8641. Four (out -16 -18 -20 -22 0 log k of 12) individual distributions do not have the power law tail at all for this threshold.
At the highest threshold (θ3), the degree distribution of the majority of networks has a power law character. The exception is one outlier. The average scaling exponent is γ3 = −1.3429 (interval [−2.1867, −0.8046].) 3
Model of the functional brain networks There are only a few papers, which attempt to model the dynamics of functional brain networks. For example Portillo and Gleiser [5] developed an adaptive complex network model, where different anatomical regions in the brain are represented by microscopic units, dynamical nodes. They start from a small random network, which grows by the addition of the new nodes with fixed number of connections. The newcomers are linked at random, but then the connections are adaptively rewired according to coherence. The state of the system is calculated at each iteration, and the evolution of nodes is given by the dynamical equation describing a set of non-linear phase oscillators. The global and local rewiring process depends on the current state. Gleiser and Spoormaker later adapted this model to model the hierarchical structure in the functional brain networks [2]. A different principle to model functional brain networks has been used by Vértes and others [6]. They proposed a model incorporating the factor of economy governing a link establishment. The topology of functional brain networks emerges from the two competitive factors: a distance penalty based on the cost of maintaining long range connections and a topological term favoring links between brain regions sharing similar input.
Model fit, participant 32, at the threshold 0.899876 -3 -4 -5 -8 -9 ) -6 k ( P g lo -7
In this paper I follow a different principle. Similar picture, as with changing the correlation threshold in our analysis, is described in Scholz et al. [7] for the noisy scale free networks. The authors started from a network with pure scale free degree distribution. Then, after fixing the number of nodes to N0, and also the initial number of edges to L0, this network is disturbed by some type of noise: i.e. random link removal, random link exchange and random link addition. The authors have studied, how the degree distribution drifts from the power law character with increasing the noise (randomness) in the network.
I observed the same pattern, namely, that the lowering of the correlation threshold is analogical to increasing the probability of addition of random links in the functional networks, which in turn causes, that the degree distribution is not power law any more. The situation can be described as follows: To utilize a view of coming nodes at each time unit, common in the growing network models, I relate time and threshold. One starts at the highest threshold θ3 (time t0 = 0), where the network is scale free having N0 nodes, L0 edges and the power law degree distribution. Then the threshold is lowered as the time flows. New nodes and edges are added to the network by both – a preferential and random linking. We suppose, that the threshold discrete and infinitesimal “jumps" can be accommodated in such a way, that only one new node and on average the same number of new edges appear per iteration. Each new node brings a1 new edges, which are linked preferentially and a edges which are linked randomly. The total amount of edges brought by one node in each time step is therefore a + a1. Simultaneously another process takes place. As the threshold decreases (time flows), some correlations between the couples of nodes already present in the network become significant. Therefore new edges are distributed randomly (b) and preferentially (b1), respectively, among the nodes being already in the network.
Thus, unlike in the model of Scholz et al [7], the network grows in the number of nodes and edges as well. Because the network at the highest threshold, which corresponds to time t0, is scale free, it is supposed, that the real correlations between voxels construct scale free structure, which is, as the threshold lowers (time grows), disturbed by the accidental correlations (links). These correlations are caused either by the real influence between the two voxels or by an accidental resemblance of the two measured signals. As we know from the theory of growing networks, scale free degree distribution is created by the preferential attachment [8].
Model fit, participant 41 at the threshold 0.819398
-4
-5
-6
-7
)(k -8
P
log -9
P(k,t + 1) = pk,k−1(t)P(k − 1,t) + (1 − pk+1,k(t))P(k,t)
(
In (
The model (
First the experimental data have been used to find the
parameters c and γ in the power law distributions at the
highest threshold θ3 (
P(k) = ckγ Both parameters c and γ are derived from the data. The power law distribution function at the highest threshold has been normalized by the constant n (based on the data) calculated from the equation n =
Z ∞ 1
P(k)dk.
and it has been checked whether the sum of all
probabilities of the initial distribution is close to one after the
normalization. The networks, for which the integral (
In the numerical simulations I first applied the model
to model the transition between the two highest
correlation thresholds, namely θ3. and θ2 of the functional brain
networks. Each model has been iterated N2–N0 times
(because at each time unit only one node appears) for the
defined set of parameters a, a1, b, b1. N0, N2 denote the
number of nodes at the initial time (threshold θ3) and at
the time t2 corresponding to the lower correlation
threshold θ2. These numbers of nodes I have from the data. In
(
Model fit for the participant 11 at the threshold 0.899876 1 2 3 4 5 6
7 log k each time step (a discrete small threshold change) a fixed number of edges is added, namely ((NL22−−LN00)) , where L2 is the number of network edges gained from the measured data at the threshold θ2 and L0 is the initial number of edges. To find the best set of parameters a, b, a1, b1 we used the hill climbing algorithm, in which the mean square error between the measured and simulated datasets has been calculated. From the best fit parameters in the current simulation fifteen new sets of parameters have been derived by slight perturbations of the currently best fit parameter set. This is a standard procedure in the hill climbing algorithm. The hill climbing algorithm has been iterated 800 times . is different, namely ((NL11−−LN00)) , where L1 is the number of edges in the functional brain network created at the lowest threshold θ1. N1, N0, L1, L0 are estimated from the data.
Because of the lack of place I present here the more
complete results for the HY group only. The best fits in
this group at the thresholds θ2, θ1 are seen at figs (
In this paper functional brain networks created at the three different correlation thresholds were analyzed. The networks have been measured in a three different groups of participants, namely the HY (healthy young), HE (healthy elderly) and AE (elderly with the Alzheimer disease) group.
Model fit participant 36 at the threshold 0.899876 Model fit, participant 4 at the threshold 0.819398 1 2 3 4 5 6 7 4 log k 1 2 3 5 6 7 8
Second, I do the same job as before to model the data at the threshold θ1. The only difference is, that the hill climbing algorithm has been iterated N1 − N0 times, where N1 is the number of nodes at the lowest threshold θ1 . Also the number of edges added in each threshold jump (time step) • There are significant intergroup differences in the degree distributions at each threshold. For example, at the highest threshold, the power law tails are less pronounced in the AE group in comparison to the HY and HE group. In one case there is no power law tail at all in the AE group. • There are also individual differences in the functional brain network degree distributions in each group at each threshold. These individual differences are most significant in the AE group. The HY group exhibits the most coherent behavior. In this group my dynamical model is most successful in fitting the data correctly.
Model A fit, participant 33 at the threshold 0.819398 4 log k
The reason of these studies was to get an insight into the in group and inter group differences in order to create an appropriate mathematical model. The overall picture was very similar to the one produced by the model of the noisy scale free network with randomly added edges Scholz et al [7].
As a second step a model of the functional brain networks has been suggested. Its detailed description is in the previous section. The model describes the dynamical processes which occur in the growing noisy, initially scale free, network. The noise in the model is due to the random distribution of a constant number of edges among a nodes being already in network. On the other hand, the preferentially distributed edges among the nodes already in network, support the scale free structure. The network also grows by the node addition, each node brings a constant number of a new edges, which are distributed either randomly or preferentially.
Due to the fact, that the original data in the HY group were the least noisy and exhibit the greatest coherence of behavior, the model gives the best results for this group. The measured data in the other groups (HE, AE) were rather noisy which influenced also preprocessing and network creation itself [3] [4]. The model gives less accurate fits in these groups, although they are qualitatively in an accordance with the data.
In conclusion, I would like to point out, that the same model accounts for the data from the HY, HE and AE individuals. This means that there might be a universal principle how the brain dynamically organizes its functional networks regardless of the age and/or onset of a neurodegenerative disease. This is a prediction arising from the model, which however needs further testing. For example the model could be enriched by taking into account a fact, that only one edge can be added at a time step to a certain node. The others should be added elsewhere. If such model will perform better, it is possible to test another one, allowing two, three... edges to be added in one time step to the same node. This is, however, left for further studies.
Regardless brain functional network, I think, that the mathematical model of growing noisy scale free network can be interesting itself as well. There might be another situations in reality to which such model can be applied.
I would like to thank prof. Benˇušková for careful reading of this text. I am also grateful to doc. Rudolf for many discussions.