ITAT 2016 Proceedings, CEUR Workshop Proceedings Vol. 1649, pp. 226–231
http://ceur-ws.org/Vol-1649, Series ISSN 1613-0073, c 2016 M. Markošová
Dynamic model of functional brain networks
Mária Markošová
Faculty of Mathematics, Physics and Informatics, Department of Applied Informatics, Comenius University, Bratislava
Abstract: Functional brain networks are networks created Functional brain networks, contrary to the structural
with a help of fMRI measurements of the in vivo brain ac- neuronal brain networks, are temporal networks. Certain
tivity [3], [4]. I elaborated a dataset, which contains func- type of the functional brain network exists only during
tional brain networks of young participants, healthy el- that time, when the brain is involved in the cognitive task
derly participants and elderly participants with diagnosed and reflects the functional cooperation of different brain
Alzheimer disease. All networks were measured at the areas. Since the smallest unit of the measured fMRI sig-
three different correlation thresholds. nal is an integrated signal of the neurons contained in one
In this paper I present a data driven mathematical model voxel, voxels are thus natural candidates for the nodes of
of functional brain networks. It is based on the threshold the functional brain network. If the two voxels function-
related shape of the degree distribution. The model is nu- ally cooperate (based on the underlying physical connec-
merically simulated and results of the simulation are com- tivity), the measured signal is highly correlated over time.
pared to the real dataset. To measure the amount of the signal correlation, the Pear-
son correlation coefficient is calculated for all the voxel
pairs:
1 Functional brain networks
r(i, j) = (1)
Functional Magnetic Resonance Imaging (fMRI) is a tech-
nique for gaining high resolution images of neural activity
< V (i,t)V ( j,t) > − < V (i,t) >< V ( j,t) >
in the brain [3], [4]). FMRI images are captured in a se- 1 1
ries of two dimensional slices, with each slice represent- (< V (i,t)2 > − < V (i,t) >2 ) 2 (< V ( j,t)2 > − < V ( j,t) >2 ) 2
ing a cross section of the brain less than 10 mm thick. A
where r(i, j) is the correlation coefficient, V (m,t) is the
single slice is comprised of a rectangular grid of discrete
measured activity in the m-th voxel at time t, and < . >
3D regions (3 × 3 × 10 mm) known as voxels (volumetric
denotes the time averages. A link between the voxel pair
pixels). A full 3D image of the brain is achieved by com-
(nodes) is established, if |r(i, j)| > θ , where θ is a pre-
bining these slices together. fMRI is an ideal technique
scribed correlation threshold. It is opted for an absolute
for deriving functional connectivity. One can ask to which
value of correlation, that is both strongly positively and
extend spatially distinct regions of the brain exhibit simi-
strongly negatively correlated voxels are included in the
lar behavior over time. By modeling this functional con-
functional network, because the functional interaction be-
nectivity as a network, we can explore the ways in which
tween neurons can be either positive (excitatory) or nega-
regions of the brain interact, and use techniques from the
tive (inhibitory). In any case, by nature, such created net-
graph theory to evaluate the topological characteristics of
works are unweighted and undirected, because correlation
these functional networks of interaction.
is a symmetric function. That means, that the node degree
In this paper I used the brain fMRI data collected by
is simply a number of the closest neighbors.
Buckner [1], (data set no. 2-2000-118W from the fMRI
Simple measures, that characterize the network in gen-
Data Center: http://www.fmridc.org). The participants
eral are averages: such as average degree, density, average
were divided into three groups: healthy young (HY) par-
shortest path, average clustering coefficient, etc. Usually,
ticipants, healthy elderly (HE) participants and elderly
these simple measures are not sufficient and one have to
participants with diagnosed Alzheimer disease (AD) (AE
rely on distributions, such as degree distribution for exam-
group). Structural and functional MRI data were acquired
ple, to acquire more detailed network properties.
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 2 Data analysis
mean age 75.1 years (SD 6.9). The AE group had 12 sub-
jects (7 females/5 males) of the mean age 77.1 years (SD In this section I complete a basic network analysis of the
5.3). There was no statistically significant difference in the measured data, already done and published by [3], [4].
mean age of the latter two groups. The networks are constructed with respect to the three
The standardized data were then used in [3], [4] to cre- different correlation thresholds, namely θ1 = 0.819398,
ate functional brain networks for each participant in all θ2 = 0.899876 and θ3 = 0.962249. McCarthy et al. [3],
three groups. These data were further elaborated by me. [4] have calculated average properties, i.e. number of
Dynamic Model of Functional Brain Networks 227
nodes, density, degree, clustering coefficient, path length, Degree distribution HY, threshold 0.899876
-1
local and global efficiency, small world index, assortativity
for all classes of functional networks. These values were -2
compared across each of the three groups, in order to find -3
gamma=-1.1407
differences related to the age and the presence of AD. I -4
concentrated my attention to the degree distributions in all
-5
three groups of participants at each of the three thresholds.
log P(k)
The reason is, that the dynamical functional brain network -6
model is based on this. -7
First, the whole brain networks of the healthy young -8
(HY) participants are described. Then I mention also the
other participant groups. At the beginning I have to state,
-9
that network is scale free if it has a power law degree dis- -10
0 1 2 3 4 5 6 7 8
tribution (4). . log k
For the lowest correlation threshold θ1 = 0.819398 one
can see in Fig. 1, that the degree distribution does not have Figure 2: HY group. Degree distribution for the functional
a power law character. Thus, the functional brain networks brain networks at the middle threshold θ2 = 0.899876.
are not scale free and the tail of the distributions is not
long enough to estimate the power law scaling exponents 0
Degree distribution HY, threshold 0.962249
correctly.
-1
-2
Degree distribution HY, threshold 0.818398
0
-3
gamma=-1.3590
-4
log P(k)
-2
-5
-6
-4
-7
log P(k)
-8
-6
-9
-8 -10
0 1 2 3 4 5 6 7
log k
-10
Figure 3: HY group. Degree distribution for the functional
0 1 2 3 4
log k
5 6 7 8
brain networks at the highest threshold θ3 = 0.962249.
Figure 1: HY group. Degree distribution for the functional
brain networks at the lowest threshold θ1 = 0.819398. For the medium threshold θ2 the degree distribution re-
veals more pronounced tail in the log–log plot and shows
For the correlation threshold θ2 = 0.899876 the degree more variability in individuals then the similar degree dis-
distribution reveals more pronounced power law tail , with tribution of the HY group. Statistical analysis of the net-
the average scaling exponent γHY 2 = −1.14 (see Fig.2). works generated for the highest threshold θ3 shows that the
The scaling exponent is calculated as an average of all degree distribution seems to be scale free, but with more
scaling exponents of all distributions for the threshold in individual differences than in the HY group. Average scal-
question. ing exponent is γHE = −1.3609 and all individual scaling
The situation changes dramatically for the highest cor- exponents are in the interval [−2.0396, −1.0500].
relation threshold θ3 = 0.962249. The functional brain We have also analyzed the functional brain networks
networks now have a well defined scale free structure, re- of the elderly people with diagnosed mild or very mild
flected in the power law degree distribution with the aver- Alzheimer disease (AE group) for all of the three thresh-
age scaling exponent γ = −1.36 (see Fig.3) with the indi- olds. In comparison to the first two groups, namely HY
vidual differences in the interval [−1.7812, −0.9346]. and HE, we have noticed greater individual differences.
The same statistical distributions as for the previous Even for the highest threshold not all of the networks are
group were analyzed for the HE group of participants. For scale free.
the lowest correlation threshold θ1 the degree distribution For the θ1 correlation threshold, the networks are not
shows similar features as in the group of the healthy young scale free, There are more pronounced power law tails of
participants, with the exception, that the individual differ- the degree distributions at the θ2 correlation threshold and
ences are more pronounced. the average scaling exponent is γ2 = −0.8641. Four (out
228 M. Markošová
Model fit, participant 29, at the threshold 0.899876 then the connections are adaptively rewired according to
-4
coherence. The state of the system is calculated at each it-
-5
eration, and the evolution of nodes is given by the dynam-
ical equation describing a set of non-linear phase oscilla-
-6
tors. 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
log P(k)
-7
brain networks [2]. A different principle to model func-
-8 tional brain networks has been used by Vértes and others
[6]. They proposed a model incorporating the factor of
-9 economy governing a link establishment. The topology
of functional brain networks emerges from the two com-
-10
0 1 2 3 4 5 6 7
petitive factors: a distance penalty based on the cost of
log k
maintaining long range connections and a topological term
favoring links between brain regions sharing similar input.
Figure 4: HY group. The best fit of the model at the
θ2 threshold. Parameters: a = 1.3957, b = 0.0013, a1 = Model fit, participant 32, at the threshold 0.899876
3.7477, b1 = 243, 1351. -3
-4
of 12) individual distributions do not have the power law -5
tail at all for this threshold.
-6
log P(k)
Model fit, participant 29, at the threshold 0.819398
-4
-7
-6
-8
-8
-9
-10
-12
log P(k)
-10
0 1 2 3 4 5 6 7 8
-14 log k
-16
Figure 6: HY group. The worst fit of the model at the
-18
θ2 threshold. Parameters: a = 1.8655, b = 0.0000, a1 =
-20 3.6148, b1 = 528, 5197.
-22
In this paper I follow a different principle. Similar pic-
0 1 2 3 4 5 6 7 8
log k
ture, as with changing the correlation threshold in our
Figure 5: HY group. The best fit of the model at the θ1 analysis, is described in Scholz et al. [7] for the noisy
threshold. Parameters: a = 12.2377, b = 1.7657, a1 = scale free networks. The authors started from a network
14.4403, b1 = 888, 3863. with pure scale free degree distribution. Then, after fix-
ing the number of nodes to N0 , and also the initial number
At the highest threshold (θ3 ), the degree distribution of of edges to L0 , this network is disturbed by some type of
the majority of networks has a power law character. The noise: i.e. random link removal, random link exchange
exception is one outlier. The average scaling exponent is and random link addition. The authors have studied, how
γ3 = −1.3429 (interval [−2.1867, −0.8046].) 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
3 Model of the functional brain networks of the correlation threshold is analogical to increasing the
probability of addition of random links in the functional
There are only a few papers, which attempt to model networks, which in turn causes, that the degree distribution
the dynamics of functional brain networks. For exam- is not power law any more. The situation can be described
ple Portillo and Gleiser [5] developed an adaptive com- as follows: To utilize a view of coming nodes at each time
plex network model, where different anatomical regions unit, common in the growing network models, I relate time
in the brain are represented by microscopic units, dynami- and threshold. One starts at the highest threshold θ3 (time
cal nodes. They start from a small random network, which t0 = 0), where the network is scale free having N0 nodes,
grows by the addition of the new nodes with fixed number L0 edges and the power law degree distribution. Then the
of connections. The newcomers are linked at random, but threshold is lowered as the time flows. New nodes and
Dynamic Model of Functional Brain Networks 229
edges are added to the network by both – a preferential In (2, 3) P(k,t) is the normalized number of nodes hav-
and random linking. We suppose, that the threshold dis- ing the degree k at the time (threshold) t. In (3) N0 , L0 de-
crete and infinitesimal “jumps" can be accommodated in note initial number of nodes and edges, a, b are the number
such a way, that only one new node and on average the of randomly added edges per iteration, where a is the num-
same number of new edges appear per iteration. Each new ber of edges fetched by a new-coming node and b is the
node brings a1 new edges, which are linked preferentially number of edges added between an older network nodes.
and a edges which are linked randomly. The total amount Similarly a1 , b1 denote the number of edges by which a
of edges brought by one node in each time step is there- new node links preferentially (a1 ) to the network and b1
fore a + a1 . Simultaneously another process takes place. is the number of edges linking older nodes preferentially.
As the threshold decreases (time flows), some correlations The transition term pk+1,k (t) describes, how the number of
between the couples of nodes already present in the net- nodes having the degree k changes due to the above men-
work become significant. Therefore new edges are dis- tioned dynamical processes. The first term of the equation
tributed randomly (b) and preferentially (b1 ), respectively, (2) is a gain term and the second one is a loss term. In the
among the nodes being already in the network. model, it is neglected what happens with the other links,
Thus, unlike in the model of Scholz et al [7], the net- the attention is payed to the fact how the link addition af-
work grows in the number of nodes and edges as well. fects the degree k (2, 3) .
Because the network at the highest threshold, which corre-
sponds to time t0 , is scale free, it is supposed, that the real
correlations between voxels construct scale free structure, 4 Results of numerical simulations
which is, as the threshold lowers (time grows), disturbed
by the accidental correlations (links). These correlations The model (3) have been simulated numerically. Each
are caused either by the real influence between the two simulation have been attempted for all functional brain
voxels or by an accidental resemblance of the two mea- networks for all the three groups of participants and com-
sured signals. As we know from the theory of growing pared to the data at the thresholds θ2 and θ1 . The best,
networks, scale free degree distribution is created by the and the worst fits for the HY group of participants and for
preferential attachment [8]. each threshold are presented here at figs 4 - 7 together with
the best fits for the HE and AE groups (figs 8 - 11). The
Model fit, participant 41 at the threshold 0.819398
networks, which were excluded, and the reasons why they
-4 were excluded, are to be explained later.
-5 First the experimental data have been used to find the
-6
parameters c and γ in the power law distributions at the
highest threshold θ3 (4). This threshold, in the threshold –
-7
time view, corresponds to the initial time t0 = 0, i.e.:
-8
log P(k)
-9 P(k) = ckγ (4)
-10 Both parameters c and γ are derived from the data. The
-11
power law distribution function at the highest threshold
has been normalized by the constant n (based on the data)
-12
calculated from the equation
-13
0 1 2 3 4 5 6 7 8 Z ∞
n= P(k)dk. (5)
log k
1
Figure 7: HY group. The worst fit of the model at the and it has been checked whether the sum of all probabili-
θ1 threshold. Parameters: a = 6.8756, b = 0.0600, a1 = ties of the initial distribution is close to one after the nor-
6.9720, b1 = 824, 4624. malization. The networks, for which the integral (5) does
not converge (in a case −1.0 < γ < 0.0 ) were excluded.
The equation describing the above mentioned dynami- Here γ is a scaling exponent of the power law degree dis-
cal processes in the model is: tribution (4) at the highest threshold.
In the numerical simulations I first applied the model
to model the transition between the two highest correla-
P(k,t + 1) = pk,k−1 (t)P(k − 1,t) + (1 − pk+1,k (t))P(k,t) tion thresholds, namely θ3 . and θ2 of the functional brain
(2) networks. Each model has been iterated N2 –N0 times (be-
In (2) the transition term pk,k−1 (t) reads: cause at each time unit only one node appears) for the de-
a + 2b (a1 + 2b1 )(k − 1) fined set of parameters a, a1 , b, b1 . N0 , N2 denote the
pk,k−1 (t) = + , (3) number of nodes at the initial time (threshold θ3 ) and at
N0 + t 2L0 + A(t)
the time t2 corresponding to the lower correlation thresh-
where A(t) = 2(a + b + a1 + b1 )t. old θ2 . These numbers of nodes I have from the data. In
230 M. Markošová
(L1 −L0 )
Model fit for the participant 11 at the threshold 0.899876 is different, namely (N , where L1 is the number of
-4 1 −N0 )
edges in the functional brain network created at the lowest
-5
threshold θ1 . N1 , N0 , L1 , L0 are estimated from the data.
Because of the lack of place I present here the more
-6 complete results for the HY group only. The best fits in
this group at the thresholds θ2 , θ1 are seen at figs (4, 5).
log P(k)
-7 The worst fits for the HY group at the threshold θ2 , θ1
are depicted at figs (6, 7). The best fits for the HE and AE
-8 groups are presented at figs 8 - 11. The comparison of the
groups is described in the discussion.
-9
-10
0 1 2 3 4 5 6 7
5 Summary and discussion
log k
In this paper functional brain networks created at the three
Figure 8: HE group. The best fit of the model at the different correlation thresholds were analyzed. The net-
θ2 threshold. Parameters: a = 2.6028, b = 0.0000, a1 = works have been measured in a three different groups of
2.5687, b1 = 313, 4685. participants, namely the HY (healthy young), HE (healthy
elderly) and AE (elderly with the Alzheimer disease)
group.
each time step (a discrete small threshold change) a fixed
(L2 −L0 )
number of edges is added, namely (N , where L2 is the Model fit participant 36 at the threshold 0.899876
2 −N0 ) -4
number of network edges gained from the measured data -4.5
at the threshold θ2 and L0 is the initial number of edges.
-5
To find the best set of parameters a, b, a1 , b1 we used
-5.5
the hill climbing algorithm, in which the mean square er-
ror between the measured and simulated datasets has been
-6
calculated. From the best fit parameters in the current sim-
-6.5
log P(k)
ulation fifteen new sets of parameters have been derived by -7
slight perturbations of the currently best fit parameter set. -7.5
This is a standard procedure in the hill climbing algorithm. -8
The hill climbing algorithm has been iterated 800 times . -8.5
-9
-9.5
Model fit, participant 4 at the threshold 0.819398 0 1 2 3 4 5 6 7
-6 log k
-7
Figure 10: AE group. The best fit of the model at the
-8
θ2 threshold. Parameters: a = 0.1516, b = 0.0000, a1 =
-9
2.0920, b1 = 228, 0264.
log P(k)
-10
As a first steps a check at which threshold the networks
-11 are scale free and how they change with the threshold
changes is performed. For this reason the degree distri-
bution for each participant at each correlation threshold
-12
-13 has been created and looked for the power law tails. It has
-14
been found found that:
0 1 2 3 4 5 6 7 8
log k
• In general, the degree distributions of the functional
brain networks changes with the threshold. At the
Figure 9: HE group. The best fit of the model at the
highest threshold θ3 the degree distributions are scale
θ1 threshold. Parameters: a = 9.3095, b = 1.7183, a1 =
free with well developed power law tail (fig. 3). As
12.3784, b1 = 3167, 9538.
the threshold decreases, the degree distributions are
changing and the power law tails are less pronounced.
Second, I do the same job as before to model the data at That means, that the network looses its scale free
the threshold θ1 . The only difference is, that the hill climb- structure (figs. 1, 2).
ing algorithm has been iterated N1 − N0 times, where N1 is
the number of nodes at the lowest threshold θ1 . Also the • There are significant intergroup differences in the de-
number of edges added in each threshold jump (time step) gree distributions at each threshold. For example, at
Dynamic Model of Functional Brain Networks 231
the highest threshold, the power law tails are less pro- fits in these groups, although they are qualitatively in an
nounced in the AE group in comparison to the HY accordance with the data.
and HE group. In one case there is no power law tail In conclusion, I would like to point out, that the same
at all in the AE group. model accounts for the data from the HY, HE and AE in-
dividuals. This means that there might be a universal prin-
• There are also individual differences in the functional ciple how the brain dynamically organizes its functional
brain network degree distributions in each group at networks regardless of the age and/or onset of a neurode-
each threshold. These individual differences are most generative disease. This is a prediction arising from the
significant in the AE group. The HY group exhibits model, which however needs further testing. For example
the most coherent behavior. In this group my dynam- the model could be enriched by taking into account a fact,
ical model is most successful in fitting the data cor- that only one edge can be added at a time step to a cer-
rectly. tain 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
Model A fit, participant 33 at the threshold 0.819398 to the same node. This is, however, left for further studies.
-5
Regardless brain functional network, I think, that the
-6
mathematical model of growing noisy scale free network
-7 can be interesting itself as well. There might be another
-8
situations in reality to which such model can be applied.
I would like to thank prof. Beňušková for careful read-
ing of this text. I am also grateful to doc. Rudolf for many
-9
log P(k)
-10 discussions.
-11
-12
References
-13
-14
[1] Buckner, R. L., Snyder, A. Z., Sanders, A. L., Raichle, M. E.,
0 1 2 3 4
log k
5 6 7 8
and Morris, J. C. Functional Brain Imaging of Young, Non-
demented, and Demented Older Adults. Journal of Cognitive
Neuroscience, 12 (Supplement 2), 24-34, (2000).
Figure 11: AE group. The best fit of the model at the
[2] Gleiser, P. M., Spoormake,r V. I., Modeling hierarchical
θ1 threshold. Parameters: a = 8.2092, b = 0.0000, a1 =
structure in functional brain network, Philos. Trans. A Math.
15.5483, b1 = 1232, 3125. Phys. Eng. Sci 368 (1933), pp 5633-44, (2010)
[3] McCarthy P., Benuskova L., Franz E. A. Functional network
The reason of these studies was to get an insight into the analysis of aging and Alzheimer disease: Results. Technical
in group and inter group differences in order to create an Report OUCS-2013-12, University of Otago, New Zealand,
appropriate mathematical model. The overall picture was (2013).
very similar to the one produced by the model of the noisy [4] McCarthy P., Benuskova L., Franz E.A. The age-related
scale free network with randomly added edges Scholz et posterior-anterior shift as revealed by voxelwise analysis of
al [7]. functional brain networks. Frontiers in Aging Neuroscience.
As a second step a model of the functional brain net- 6:301. doi: 10.3389/fnagi.2014.00301, (2014)
works has been suggested. Its detailed description is in [5] Portillo, I. J. G., Gleise,r P. M, An adaptive complex network
the previous section. The model describes the dynamical model for brain functional networks, PLoS One 4 (9), pp 1-
processes which occur in the growing noisy, initially scale 8, (2009) .
free, network. The noise in the model is due to the ran- [6] Vértes, P. E., Alexander - Bloch, A. F., Gogtay N., Giedd, J.
dom distribution of a constant number of edges among a N., Rapoport, J. L., Bullmore, E. T., Simple models of hu-
nodes being already in network. On the other hand, the man brain functional networks, Proc. Natl. Acad. Sci. USA
preferentially distributed edges among the nodes already 109 (15), pp. 5868-73, (2012)
in network, support the scale free structure. The network [7] Scholz, J., Dejori, M., Stetter, M., Greiner, M., Noisy scale
also grows by the node addition, each node brings a con- free networks, Physica A 350 (2-4), pp 622 -642, (2005)
stant number of a new edges, which are distributed either [8] A. L. Barabási, R. Albert, mergence of scaling in random
randomly or preferentially. networks, Science 286, 3616 (1999)
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 net-
work creation itself [3] [4]. The model gives less accurate