Dynamic fMRI brain connectivity: A study of the brain’s large-scale network dynamics

(1)

Master’s Thesis

Dynamic fMRI brain connectivity:

A study of the brain’s large-scale network dynamics

Author: Supervisors:

Per Brantefors William Hedley Thompson

pebr0018@student.umu.se Peter Fransson Examiner:

Anders Eklund

Karolinska Institutet

Department of Clinical Neurosience Master’s Thesis Engineering Physics

Academic Year 2015/2016, January 18, 2016

(2)

Abstract

Approximately 20% of the body’s energy consumption is ongoingly consumed by the brain, where the main part is due to the neural activity, which is only increased slightly when doing a demanding task. This ongoingly neural activity are studied with the so called resting-state fMRI, which mean that the neural activity in the brain is measured for participants with no specific task. These studies have been useful to understand the neural function and how the neural networks are constructed and cooperate. This have also been helpful in several clinical research, for example have differences been identified between bipolar disorder and major depressive disorder [1]. Recent research has focused on temporal properties of the ongoing activity and it is well known that neural activity occurs in bursts. In this study, resting-state fMRI data and temporal graph theory is used to develop a point based method (PBM) to quantify these bursts at a nodal level. By doing this, the bursty pattern can be further investigated and the nodes showing the most bursty pattern (i.e hubs) can be identified. The method developed shows a robustness regarding several different aspects. In the method is two different variance threshold algorithms suggested. One local variance threshold (LVT) based on the individual variance of the edge time-series and one global variance threshold (GVT) based on the variance of all edges time-series, where the GVT shows the highest robustness. However, the choice of threshold needs to be adapted for the aims of the current study. Finally, this method ends up in a new measure to quantify this bursty pattern named bursty centrality. The derived temporal graph theoretical measure was correlated with traditional static graph properties used in resting state and showed a low but significant correlation. By applying this method on resting-state fMRI data for 32 young adults was it possible to identify regions of the brain that showed the most dynamic properties, these regions differed between the two thresholding algorithms.

(3)

Sammanfattning

Ungefär 20% av kroppens energikonsumtion konsumeras ständigt av hjärnan, där huvudelen av enegin går åt till den neurala aktiviteten, som bara ökar nå- got vid en pågående arbetesuppgift. Denna pågående neurala aktiviteten stud- eras med fMRI i vilande tillstånd, vilket innebär att den neurala aktiviteten i hjärnan mäts när testdeltagaren inte har någon specifik uppgift. Dessa studier har varit viktiga för att få förståelse för de neurala funktionerna och hur de neurala nätverken är uppbyggda och sammarbetar. Detta har också varit till nytta i flertalet kliniska projekt, till exempel har skillnader identifierats mellan bipolära personer och personer med allvarlig depression [1]. Senare forskning har fokuserat på de temporala egenskaperna av den pågående neurala aktiviteten och det är nu allmänt känt att den neurala aktiviteten sker i skurar. I den här studien används data från fMRI i vilande tillstånd och temporal grafteori att för att utveckla en punkt baserad metod (PBM) för att kvantifiera dessa skurar på nodnivå. Genom denna metod kan man ytterligare undersöka detta mönster och de noder som uppvisar den högsta graden av skurar (s.k hubar) kan identifieras.

Metoden som har utvecklats uppvisar en robusthet i flera olika avseenden. I metoden presenteras även två olika algoritmer för varianströskling. En lokal var- ianströskel (LVT) som baseras på den individuella variansen för varje kopplings tidsserie och en global varianströskel (GVT) som utgår från den globala variansen för alla kopplingars tidsserier, där GVT uppvisar den bästa robustheten.

Men valet av trösklings algoritm måste anpassas för målet man har med studien.

Slutligen presenteras ett nytt mått för att kvantifiera detta mönster på nodnivå, bursty centrality. Det utvecklade temporala grafiteori måttet korrelerades sedan med traditionella statiska grafteori mått som ofta används för fMRI i vilande tillstånd och visade en svag men signifikant korrelation. Genom att sedan an- vända denna metod på data från fMRI i vilande tillstånd från 32 unga vuxna var det möjligt att identifiera områden i hjärnan som visade de största dynamiska egenskaperna, dessa områden varierade för de olika trökslings algoritmerna.

(4)

(5)

Commonly Used Acronyms

MRI Magnetic Resonance Imaging

fMRI functional Magnetic Resonance Imaging

TR Repetition Time; indicating the temporal resolution of MRI images

Hb Oxygenated hemoglobin dHb Deoxygenated hemoglobin BOLD Blood-Oxygen-Level-Dependent

MNI Atlas of brain coordinates produced by Montreal Neurolog- ical Institute

AAL Automated Anatomical Labeling

BA Brodmann area

PBM Point based method ROI Region of interest Power Power 164 atlas HO Harvard-Oxford atlas

PCA Principle component analysis DMN Default mode network SMN Somato-motor network Vis Visual network

FPN Frontal parietal network SAN Salience network

CON Cingulo opercular network AU Auditory network

Sub Subcortical network VAN Ventral attention network DAN Dorsal attention network ETPS Edge time-point scaling GTPS Global time-point scaling LVT Local variance threshold GVT Global variance threshold ISI Inter spike interval

min_L Minimum length of a burst

(6)

Dynamic fMRI brain connectivity: Contents

1 Introduction

1.1 Background

There are several ways to investigate the neural activity in the brain, for example electroencephalography (EEG), magnetoencephalography (MEG) and the one used in this project, functional magnetic resonance imaging (fMRI).

In the beginning, the majority of the fMRI studies in neuroscience have been focusing on the brain response to a certain task or a stimuli [2]. These studies have linked certain regions of the brain to specific tasks or a specific stimuli.

Based on responses to tasks and stimuli, the neural structure have been divided into several functional networks. It has also been noted that there is some intrinsic fluctuating neural activity in the brain even during resting state conditions (i.e when no task are performed and no stimuli is applied), which is the reason for the resting-state fMRI studies. These resting-state studies have identified a region of brain activities active during resting conditions, which has been named the default mode network (DMN). Before Biswal in 1995 first studied these spontaneous fluctuations in neural activity [3] were these considered to be noise, but due to several reasons, these fluctuations is important to study further. Firstly, the brain accounts for approximately 20% of the total energy consumption during resting state but it weighs only about 2% of the body weight. Most of the energy consumption in the brain is due to neural activity [2].

The increase in metabolism during a task is usually less then 5% compared to the resting state consumption. This increase is very small because the activity is increased in a very small part of the brain. To get a better knowledge of how the brain operates, the spontaneous neural activity has to be studied, which is what is isolated during rest. Studies have also shown that these spontaneous activities are not random, one example of this is the higher activation in the DMN during resting state, another is the fact that left somatomotor cortex is specifically correlated with the right somatomotor cortex during resting state [3]. The neural fluctuation during resting-state condition is therefore important to study further, which was done for 10 resting-state networks in this study.

In earlier studies have a bursty pattern in the neural activity been suggested [4]. It was also suggested that the bursty pattern are stronger between networks while the connectivity within a network displays both a bursty and a periodic pattern of connectivity [4]. However, these studies only shows that there is a bursty pattern in the neural activity and not which nodes or networks that are the most bursty ones. Now, the temporal properties of resting-state connectivity is being investigated, for example the fact that the connectivity increases in a bursty pattern in large scale connectivity. To do this, temporal graph theory have been applied were the brain is seen as a mathematical system consisting of several items called nodes and connection between them called edges and a set of nodes and edges is called network. This will hopefully identify certain regions (i.e nodes) that are more important for the dissemination of information in the brain. Further, also a simultaneously activation in different nodes indicate that the information is sent between these nodes. This can indicate whether there is a specific node that is important for sending information within the network respectively between the networks. Once a methods to quantify the dynamic nature of the spontaneous neural activity are in place, they can be applied in, for example, clinical settings where traditionally resting state activity has been use [5]. However, to do this, the first thing needed to be done is to develop a

(9)

Dynamic fMRI brain connectivity: 1 Introduction

method and a measure for this quantification on a smaller level (i.e node level) in the brain.

1.2 Aims

The main goal of this project is to in a graph theoretical way study the bursty pattern in the resting-state fMRI data shown by Thompson and Fransson in the article Bursty and persistent properties of a large-scale brain networks revealed with a point-based method for dynamic functional connectivity [4]. The method used in the said article will form the foundation of this work but will be extended in several ways. In the end, the main goal is to quantify the burstiness on a nodal level by introducing a new graph theoretical measure. The purpose of this measure is to quantify which nodes and which networks that shows the most bursty pattern and thus have an important role in the dissemination of information in the brain.

1.3 Limitations

All the analysis in this project is limited to a graph theoretical approach and the development of the method will be performed only in Matlab and Matlab-based programs. The data used to develop the method are limited to resting-state fMRI eyes open data of 56 healthy adults found in the Oulu A dataset, but also the Oulu B dataset including 47 healthy adults will be used in the end to confirm the results.

(10)

Dynamic fMRI brain connectivity: 2 Theory

2 Theory

2.1 Imaging

2.1.1 Magnetic Resonance Imaging

Functional Magnetic Resonance Imaging (fMRI) is a powerful tool to measure the neural activity in the brain. It is based on Magnetic Resonance Imaging (MRI) which uses the properties of the hydrogen nucleus to create a image of an object. The body consists of a large amount of water molecules and each water molecule contain two hydrogen nuclei. Since the hydrogen nucleus consists of a single proton which have an intrinsic spin, a magnetic field is produced, also called a magnetic moment. Due to this, all the hydrogen atoms will align to the strong external magnetic field (B₀) applied in the MRI-scanner [6] according to the Bloch equation

dM(t)

dt = γM(t) × B(t) (1)

where M is the magnetization and γ is the gyromagnetic ratio. The gyromagnetic ratio for protons is equal to 2.7 × 10⁸ rad/sT or 42.57 M HzT⁻¹ [7]. In fact, the protons aligned with the magnetic field will also be rotating slightly around the external magnetic field. This is called precessing, and the pressesion rate is termed Larmor frequency, ωL, see figure 1. The Larmor frequency is related to the magnetic field according to the Larmor equation below

ω_L= γB₀ (2)

When a radio frequency pulse at the Larmor frequency is applied, the atoms will flip into a higher energy state and the magnetization of the hydrogen nuclei will tilt to the transverse plan orthogonal to B₀. The magnetization will as a reaction go back to the equilibrium state along B0. Thus, the absorbed energy will be emitted as radio frequency signal when flipping back to the low energy state [6]. The time it takes for the magnetization to align to B0 after the excitation is tissue specific and called spin-lattice relaxation, T1. Initially all the protons precesses will be in phase. But due to interactions between the hydrogen atoms the protons will begin to dephase, this leads to a decrease in magnetization over time and is known as the spin-spin relaxation, T2. But due to inhomogeneities in the magnetic field and the magnetic susceptibility an additional dephasing called T₂^∗ will occur. The magnitude of the radio frequency waves emitted during the relaxation is proportional to the transverse component of the net magnetization and is detected by the radio frequency coil in the scanner.

To localize the signal, three magnetic gradients are used. By applying a gradient field, a linear variation in the magnetization is created which as a result will produce variations in the resonance frequency in the direction of the gradient field. This means that the radio waves with a certain frequency only will fulfill the resonance condition in equation 2 at a certain plane perpendicular to the gradient field, this is called slice selection [6]. To localize the signal in the slice the so called phase encoding gradient and the frequency encoding gradients are applied. The phase encoding gradient rendering different phases of the magnetization and the frequency encoding gradient rendering different frequencies in the slice [8]. By combining the phase and frequency at a specific time each spatial location in the slice can be localized.

(11)

By changing the duration and onset times for the gradients, the radio frequency pulse and the time-echo (TE), which is the time between the magnetization is flipped until the signal is measured, different sequences can be obtained.

By using different sequences is it possible to obtain different image contrasts.

Further more, the contrast in the image also depends on the tissue specific relaxation. To imaging the blood oxygen level dependent (BOLD) signal, which is of interest in this case, a so called echo planar imaging (EPI) sequence is used.

The EPI is a T₂^∗ weighted sequence.

Figure 1: Precessing of the hydrogen proton at the larmor frequency (ωL) around the main magnetic field (B0). (Image source: https : //en.wikipedia.org/wiki/Larmor_precession)

2.1.2 Functional Magnetic Resonance Imaging

There is about 100 billion neurons in the brain [9], so how do we measure the activity of those? There are several techniques to measure the activity non- invasively in humans, for example electroencephalography (EEG) and magnetoencephalography (MEG) which measure the electrical and the magnetic signals caused by the neural activity using several sensors applied to or near the head [9]. However, these techniques have a poor spatial resolution and also other limitations. fMRI on the other hand is an indirect measurement technique which measure the neurovascular changes that are correlated with the neural activity.

The neural activity require energy to work, in fact the brain it self consumes around 20% of the total energy consumption in the body [2]. It is the vascular system that supply the brain with the required energy through the delivery of oxygen and glucose in the blood. The fact is that the hemoglobin has differences in magnetic properties depending on whether it is bounded to oxygen or not. Oxygenated hemoglobin (Hb) is diamagnetic which means that it has no unpaired electrons and hence no magnetic moment. Deoxygenated hemoglobin (dHb) on the other hand has both unpaired electrons and a significant magnetic moment and is hence paramagnetic [9]. The difference in magnetic susceptibility between oxygenated and deoxygenated hemoglobin can be around 20% greater

(12)

for the deoxygenated hemoglobin. This is the physological property which is used in the fMRI scanner [9]. The paramagnetic hemoglobin will affect the magnetic field, also nearby protons will be affected and as a result precess at different frequencies. This will result in a shorter T₂^∗ which leads to the oxygenated blood will have a stronger MR signal for a T₂^∗ sensitive sequence than the deoxygenated blood will have. This signal is called the blood-oxygen-level- dependent (BOLD) contrast.

2.1.3 Spatial resolution

The spatial resolution of fMRI data is often stated as voxel size which is determined by the three scanner parameters: field of view, matrix size and slice thickness. When studying the entire brain using fMRI a voxel size of around 4 × 4 × 4 mm is common, while a quite smaller one, approximately 1 × 1 × 1 mm, is used when studying a single region of the brain. For anatomical images a voxel size of approximately 1 × 1 × 1 mm is used [9]. To make it easier to visualize the results is it common to use a anatomical scan to display the functional data upon. So why do one often use such a large voxel size for the functional data? There are mainly two reasons. First, a smaller voxel size will decrease the signal-to-noise-ratio of the BOLD signal. This is because the total amount of deoxygenated hemoglobin within a voxel affect the BOLD signal. This means that when reducing the voxel size by half also the BOLD signal changes will be reduced by half, which results in a smaller signal which leads to a smaller signal-to-noise ratio [9]. The second reason for not using a smaller voxel size is the scanning time, which leads to a reduced temporal resolution. For example when decreasing the slice thickness by half, twice as many slices is needed which of course will double the acquisition time. Also the slice resolution affect the scanning time. A typical scanner, for example, could acquire 16 slides per second at a resolution of 4 × 4 × 4 mm and only 1 slice per second at 1 × 1 × 1 mm resolution.

2.1.4 Temporal resolution

The temporal resolution for fMRI is usually expressed as repetition time, or TR.

The neurons is active less then a second, which would require a high temporal resolution. But since fMRI measure slower changes in the vascular system i.e.

the BOLD signal, which is active for approximately 10 seconds following some event, the temporal resolution can be fairly low. However, increasing the temporal resolution will of course improve the estimates of the vascular changes. So the temporal resolution is not only dependent on the TR but also by the limitations in the vascular system. However, a too short TR can also cause some problems. For example is the spatial resolution decreased with a shorter TR, which we touched upon in section 2.1.3. It is also possible to acquire stronger MR signals when using a larger TR, this is because when the RF pulse tilt the magnetization it must have enough time to reach its equilibrium position before the next RF pulse. So if a too short TR is used the tilt angle have to be reduced, and since the MR signal is proportional to the projection in the transverse plane also the MR signal will decrease. In fMRI studies a TR around 1.5-3 seconds is common [9].

(13)

2.2 Neuroanatomy and neural networks

2.2.1 Neurons

The central nervous system (CNS) consists of the spinal cord and the brain.

Further, the CNS is composed of several cellular elements including the neurons which are the principal information processing cells. The neurons consists of cell bodies, axons and dendrites. The axon is the connection to other neurons and the dendrites is the connection to other cells. The brain can be divided into three different structures; gray matter, which primarily consists of cell bodies, white matter, which primarily consists of axons and the cerebrospinal fluid (CSF), which is a colorless liquid [9]. Since the BOLD signal is derived from the gray matter is this the interesting part when studying the neural activity using fMRI.

Anatomists have divided the brain into several different regions in many different ways, two examples are the Automated Anatomical Labeling (AAL) and the Brodmann areas (BA). The Brodmann areas was developed by Korbinian Brodmann in 1909 who divided the brain in 47 different regions and have since that been developed to 52 regions. The BA is based on the composition and structure of the neurons in the cerebral cortex, which means that it is a structural partition. AAL on the other hand is a anatomical partition divided into 90 different regions based on the structure of the brain [10]. An example of the different partition, both for AAL and BA can be seen in figure 2.

Figure 2: Sectional image of the AAL and Brodmann parcellation A) AAL partition in the coronal plane. B) AAL partition in the transverse plane. C) AAL partition in the sagittal plane. D) Brodmann partition in the coronal plane. E) Brodmann partition in the transverse plane. F) Brodmann partition in the sagittal plane.

2.2.2 Neural networks

The brains anatomical areas are connected in networks. Correlated brain activity on these anatomical networks create what is called functional connectivity.

These functional networks are often isolated during periods of rest and called

(14)

resting state networks (RSNs). In this study, these networks are derived using a template (Power 264 atlas, see section 3.6), into 10 RSNs. Since the purpose of this study together with many other studies is to mapping the brain’s functional network, there is still some uncertainties exactly what each of these neural networks do but below is a short summary of each of the 10 resting-state networks used in this study. In figure 3 an illustration of the networks is given.

Somato-Motor (SM)

The first identified networks and is mainly linked to sensory perception and movement of the body.

Cingulo-Opercular (CO)

Have been hard to functionally characterize but it seems to bee related to the maintenance activity that spans over the whole task [11].

Auditory (AU)

Linked to hearing and processing sound input.

Default mode network (DMN)

Connected to self and social cognition and is the part of the brain that are hypothesized to be inactive during a task and are more active during resting state conditions [12].

Visual (VIS)

Linked to process visual information.

Frontal parietal (FPN)

Supports phasic components of attentional control (adaption after errors, task-switching).

Salience (SA)

Supports introspective awareness and autonomic processes, for example breathing, blood circulation etc.

Subcortical (SUB)

A collection of the areas which not is located in cortex. These areas is important for biological functions such as hunger and temperature regu- lation. It also have an important role for relaying the incoming signals to the other parts of the brain.

Ventral Attention (VA)

Related to unvoluntary attention and consciousness. This area supports the detecting of unexpected or unattended stimuli.

Dorsal Attention (DA)

Related to voluntary attention and consciousness and supports selective attention in the visual and spatial domains [13].

Undefined (U)

Contains nodes that could not be included in any of the networks or pos- sibly belong to multiple networks and therefore were hard to classify.

(15)

Figure 3: The resting-state network partitions used in this study.

2.3 Graph theory in neuroscience

To study resting-state fMRI many different methods can be used. One of the most powerful, and the one used in this project is the graph theoretical approach.

In the graph theoretical approach a mathematical system consisting of items and pairwise relationship between the items describes a complex system. The items are called nodes, the relationships between them are called edges and a set of nodes and their edges are called a network [14]. From these is often a so called connectivity matrix, Ai,j (or adjacency matrix) created, where i and j is nodes and A representing the connectivity between the nodes (i.e edges), see figure 4.

By using this approach is it possible to identify and quantify substructures and hierarchy within a network, identify hubs and critical nodes and see how the information in the network is distributed [14]. Another main advantage of using graph theory is the simplicity to model the systems in different levels such as the entire graph, subgraphs (i.e. networks) and individual nodes. An extension of graph theory is quantifying the changes in connectivity during time, temporal graph theory [15].

Graph theory is useful to consider the connections of multiple brain areas. And also use different graph theoretical measures, for example betweenness centrality, C_B. Betweenness centrality is a measure of the fraction of shortest paths in a network that contain a specific node. A high betweenness centrality value thus means that the current node participate in many shortest paths [15], and thus indicates that the node is important for information dissemination. Betweenness

(16)

centrality can also be described with the following formula.

CB(i) =Σi6=j6=kvi(j, k)

Σi6=j6=kv(j, k) (3)

Where v(j, k) is the number of paths between j and k and vi(j, k) is the number of shortest paths between k and j that goes through i.

In the work by Thompson and Fransson [4] the temporal graph theoretical measure burstiness was studied for resting-state fMRI data. Burstiness measure if the interaction between two nodes appears in bursts, which means that the interactions between two nodes are characterized by short periods of time with much interaction followed by longer periods of time without any interactions.

To see whether the information is sent in bursts the measure in equation 4 was used [4].

B = σ_τ− µ_τ στ+ µτ

(4) where σ_τ are the standard deviation and µ_τ is the mean of the inter contact times. From equation 4 one can see that B will range between -1 and 1. If B > 0 a bursty pattern is suggested and if B < 0 a more periodic pattern is suggested. The study by Thompson and Fransson [4] showed that the connectivity profiles contained a mixture of bursty and periodic pattern. However, the bursty patterns was dominating mainly between the resting state networks [4].

This thesis is mainly based on the work by Thompson and Fransson [4] and will further study this bursty pattern on a nodal level. Some parts of the method used in this thesis is also introduced in said article.

Figure 4: Connectivity matrix, where the nodes are represented on the axes and the relationship between the nodes (i.e. the connectivity) is represented at each pixel. In this case is the nodes sorted due to their network affiliation which makes it possible to see the higher connectivity within the networks as squares along the diagonal.

(17)

Dynamic fMRI brain connectivity: 3 Method

3 Method

3.1 Data

The analysis was performed on resting-state fMRI data included in the Oulu A data set taken from the 1000 functional connectomes [16]. The data was collected during resting-state with eyes open using a fixation cross for 56 healthy subjects at the age ranging 20 − 22 years old. There were 23 males and 33 females. The images included the cerebellum which sometimes led to the top parts of the brain being missing which led to some of the subject had to be rejected after the preprocessing (see section 3.2) due to incomplete data. A total of 24 subjects were rejected yielding that 32 subjects remained in the study.

The age of the remaining subjects ranged from 21 − 22 years old including 8 males and 24 females.

To verify the results was also the Oulu B data set used, also taken from the 1000 functional connectomes [16]. This consists of 47 healthy subjects at the age ranging 20 − 23 years old. There was 14 males and 33 females. Due to the same reason as for Oulu A was 18 subjects rejected yielding that 29 subjects remained. The age of the remaining subject ranged from 20 − 22 old, there was 8 males and 21 females.

The fMRI data was in both cases collected on a 1.5T scanner with an repetition time (TR) of 1800 ms and 28 oblique axial slices at 245 time-points. This data was selected since it was found to have little sleep artifacts which can be a problem with resting state data [17].

3.2 Preprocessing - SPM12

To assemble the functional and anatomical data and to make the different subjects comparable, preprocessing was necessary. This was done by the Matlab based program SPM12. The preprocessing was done in five steps, head motion correction, coregistration, segmentation, normalisation and smoothing. Below will these steps be described in brief, for more information about SPM12 see SPM12 Manual [18].

3.2.1 Head motion correction

The most important step in the preprocessing is to correct for head motions during the scan. Since the voxels are defined by their position in the scanner and not their position in the brain, head movements can result in some huge errors if they are not corrected. For example, if the subject moves his or her head inside the scanner the scanner will record the activation in a specific spatial locations in different voxels at different time-points.

Since the head movement in a fMRI scan often are quite small, the change in shape or size are also very small and can therefore be ignored. Because of this SPM12 consider the head to be a rigid body. The motions of this rigid body can be described by six variables, z which is defined along the magnetic field, thus runs from the feet to the top of the head of the subject. The x-axis

(18)

runs through the subjects ears, from left to right, and the y-axis runs from the back of the head and exits the forehead. Yaw is the rotation around the z-axis, pitch the rotation around the x-axis and roll is the rotation around the y-axis.

To do these head motion corrections SPM12 uses the BOLD response of one image, selected by the user, as reference and then correct all the other images at each time-points as good as possible at each coordinate. This assumes that the BOLD signal is exactly the same for all time-points which of course not is the case. However we expect that the activation in the most voxels will be roughly the same for the functional data through all the time-points which is used to do the correction as good as possible [19].

3.2.2 Coregistration

During the functional scan the whole brain is scanned once every 1.8 seconds to obtain a useful temporal resolution. Due to this short scanning time the spatial resolution become poor, in this case, the functional data have a voxel size of 4 × 4 × 4.4 mm. In the structural scan may it take several minutes to scan the whole brain once with a voxel size of 1 × 1 × 1 mm. Since this implies that the spatial resolution is lower for the functional images compared to the structural images is it needed to align the functional and the structural images, this is called coregistration. The coregistration also makes it possible to improve the spatial localization of the functional data with the use of the spatial resolution of the structural image. Since the functional and the structural images have different voxel sizes and since the contrast often are different in the two images the most common approaches like sum of squared difference will not work [19].

Instead a frequency histogram of the intensity values in each image is created, then the intensity values in every voxel is replaced by the bin number from the histogram which corresponds to that voxels intensity value. To finally compare the association between the two images the measure mutual information from information theory is used and the rigid body movement that maximize the mutual information between the histograms are found.

3.2.3 Segmentation

Since the data contains the information of the whole brain but BOLD signal is derived only from the gray matter are the images segmented. The segmentation divides the images in mainly four parts, these parts are gray matter, white matter, cerebrospinal fluid and bone structure. The segmentation algorithm in SPM12 is based on the algorithm described by John Ashburner and Karl J.

Friston in the paper Unified segmentaion [20] with some minor changes [18].

3.2.4 Normalization

Since there are large differences in shape and size between individual brains, there are also large differences in every brain region. This makes it hard to know exactly were in the neuroanatomical brain structure a certain activation occur.

To mitigate these problems SPM12 register the structural images for all subject to a standard brain where the coordinates for all brain regions already are published in an atlas. In this case, the MNI (Montreal Neurological Institute)

(19)

atlas have been used. The MNI atlas is an average of 152 different brains [19].

To do this the first step is to align the image so that it will be in the same location, orientation and size as the reference, this is done by a linear affine transformation [19]. The affine transformation include translation, rotation, scaling and shearing transformations. The main difference remaining now is due to local differences between the image and the reference. To align these differences a non-linear transformation is used to perform local changes like shrinking, stretching, spatial movements etc. However, due to these changes all the voxels will no longer be the same size and not rectangular and therefore it is necessary to redefine the voxels. To do this, interpolation is required to estimate the intensity value in the redefined voxels. For further information about the normalization algorithm see the SPM12 Manual [18] or the article Unified segmentation [20] which presenting the algorithm SPM12 used as a base for their algorithm. Due to the normalization is it now possible to report and interpret spatial locations in a consistent manner across subjects. Results can be compared across studies, averaged across subjects and generalized to larger populations. But there are of course also some disadvantages, for example is the spatial resolution reduced and errors may be introduced in the data.

3.2.5 Smoothing

Since the fMRI data is noisy and the change in the BOLD signal is small the signal-to-noise ratio is low. To decrease the noise and thereby increase the signal-to-noise ratio, the data is spatially smoothed. In this case the images is smoothed with a 8 × 8 × 8 mm gaussian filter kernel with the width defined by the full width at half maximum (FWHM). Beside the increase in signal-to- noise ratio, there is also other advantages with smoothing the fMRI data. For example, the data become more normally distributed when smoothed which is an advantage since some of the dominating models for fMRI data analysis assume normal distributed noise [19].

3.3 Postprocessing - CONN

After the preprocessing some extra processing was needed, this was done in the Matlab based program CONN [21]. CONN is a software designed for analysis, display and computation of functional connectivity for fMRI data. In this project was CONN mainly used for some minor processing steps. First, contributions from confounders was regressed out, for example was the white matter, cerebrospinal fluid and some movement regressors removed. Further, the signal were filtered with a bandpass filter with the cut-off frequencies 0.01 - 0.1 Hz.

Thus, contributions to the BOLD signal originating from for example heartbeat, breathing etc was reduced. The data was also detrended and finally scrubbed, this means that all of the 6 motion parameters are converted into a time-series and then is the time-series thresholded with a threshold of 0.5mm which corresponding to the 95th percentiles in a normative sample. This additional motion correction is called scrubbing and is done due to micromovements that can effect fMRI resting-state connectivity [22].

(20)

3.4 MRIcron

As was said in section 2.2.1 can the brain be parcellated into several subparts.

In this project two of the most common parcellations was used, the Automated Anatomical Labeling (AAL) and Brodmann area (BA). MRIcron is a visualization and region drawing software for neuroimaging developed by Chris Rorden [23] which uses the MNI coordinates to look up which of the AAL or Brodmann areas a specific node belongs to. This can be used to identify specific brain regions (e.g. which brain region a node from Power 264 atlas belongs to, see section 3.6). It also makes the results easier to comprehend and to compare with other studies which may have used other sets of nodes or even not a point based method which is used in this study and will be further described in the section 3.6.

3.5 BrainNet Viewer

BrainNet Viewer [24] is a simple Matlab-based program for visualization of structural and functional brain networks and is suitable when using an graph theoretical approach. In BrainNet Viewer, a template of a brain is loaded as a surface and further the brain connectivity data is loaded. The data can be loaded as a volume, nodes or as edges. There is also possible to do some plotting settings in BrainNet Viewer such as thresholding, color coding, size coding etc.

This makes it possible to visualize the results in a distinctly and easy way.

3.6 Power 264 atlas

After the processing steps, the nodes was created. One standard method is simply define each voxel as a node. But since this project investigate functional areas on a macroscopic level a point based method was used, this means that several region of interests (ROI) is extracted. The set of ROIs extracted is the set Power et al first defined in the article Functional Network Organization of the Human Brain [14]. This gives 264 ROIs covering the cerebral cortex, subcortical structures and cerebellum. Each ROI was in this case a sphere with a 10 mm radius which was defined as one node. Since the ROI cover several voxels, the mean value of the BOLD signals for all these voxels was used. This will of course reduce the spatial resolution. The reason to extract the ROIs using Power 264 atlas is because it also divided the ROIs into the networks described in section 2.2.2 based on the connectivity between them [14]. All the ROIs are shown in figure 5 where the color indicates their network affiliation, these can also be seen in figure 3 in section 3. Out of the total 264 ROIs is 58 in the DMN, 35 in SM, 31 in Vis, 24 in FP, 18 in Sa, 14 in CO, 13 in Au, 13 in Sub, 11 in DA, 9 in VA and 38 nodes are undefined with respect to resting-state network.

(21)

Figure 5: Parcellation of the brain using the Power 264 nodes where the color shows there network affiliation, the white nodes is the ones with an undefined network affiliation. (Image source: Bursty and persistent properties of a large- scale brain networks revealed with a point-based method for dynamic functional connectivity [4])

3.7 Harvard-Oxford atlas

To check how the results depends of the chosen Power 264 atlas described above, also another atlas was used during some parts of the project. This atlas is called the Harvard-Oxford atlas (HO-atlas). The HO-atlas cover the whole cortex and also some subcortical regions and consists of 57 brain areas but splitting many of the structures into left and right hemisphere ended up in 106 different regions. This parcellation was taken from CONN toolbox. To make it easier to evaluate the results, these region were also divided into networks by the use of the code infomap developed by Roswall and Bergström and available at http://www.mapequation.org/code.html [25]. This resulted in 17 different sub- networks where many of them contained only a few number of regions. All of the networks containing less then five regions were put together into one network, this resulted in a total of 7 networks which can be seen in figure 6 below.

Figure 6: The network parcellation for the Harvard Oxford atlas. This consists of 7 networks, showed by the color

3.8 Principal Component Analysis

Principal component analysis (PCA) is a multivariate statistical technique and a common technique for data reduction. PCA finds a linear transform which

(22)

represent as much as possible of the variance in the data in fewer dimensions.

The bases of the transformed data is the eigenvectors and the eigenvalues represent the variance [19]. When the variance is represented in a few dimensions, the new data in the PCA dimension is sorted starting with the one containing the highest variance. By then keeping the most prominent PCA dimensions, i.e the PCA dimensions with the highest variance, the data can be reduced but still keep much of the important information. PCA can also be used to reduce noise since the noise often are characterized by a low variance. This is done by only retain very few of the PCA dimensions with the highest variance but it also result in dynamics to be lost.

In this case data reduction was the purpose of the PCA since k-means clustering (see section 3.9) will be more efficient with a lower number of dimensions in the data. All of the 264 ROIs time-series was transformed into PCA-space and described along 264 PCA dimensions sorted due to how much of the variance each PCA dimension explains of the data. In this case, 90% of the variance, corresponding to 104 out of the 264 PCA dimensions was kept. These 104 PCA dimensions was later used in the k-means clustering in section 3.9. From the PCA also a weight called loadings is obtained, this weight describes the original data in terms of the PCA-dimensions, which is described in equation 5.

y = xl (5)

Where y is the original data x is the data in PCA-space and l is the loadings.

This will be further used in section 3.11.2.

3.9 K-means clustering

Each ROI is now represented along 104 PCA dimensions which is used in the k- means clustering. The aim of this is to divide the BOLD fMRI signal time-series into several clusters based on their co-activation pattern in their spatial dimensions regardless of their temporal location. An example for a two node system of this can be seen in figure 8A. K-mean is a iterative clustering algorithm commonly used in neuroscience [26], which divide the data into k number of cluster, where k is decided by the user. The method first places k number of centroids at random locations. Then all the data points, in this case all the time-points in each of the 104 PCA-dimension, finds the nearest centroid and assign the time-point to that specific cluster. A new location for all the centroids is then calculated by taking the mean position of all the data-points within the cluster.

Since the centroid now is in a new position new time-points will be assigned to that cluster. This is then repeated until all the clusters are unchanged. To get an even better result, independent of the initial position of the centroids, all this is then replicated several times, every time with a new random initial position for the centroids. Finally the result with the lowest sum of point-to-centroid distances is selected as the final result. There are also several other ways to calculate the new position of the centroids, the one described above is called the squared Euclidian distance and is considered to be the standard method.

The k-mean clustering was performed for k-values between 1 and 30. The clustering was performed up to 1000 iterations for each k-value to attain convergence and each k-value was repeated 20 times to reduce the influence of the initial condition and to obtain the best possible result. The centroid position from the different data-points was calculated using the squared Euclidian distance.

(23)

Finally to know which k-value to use equation 6 was used.

s(i) = b(i) − a(i)

a(i) + b(i) (6)

Here a(i) is the distance between each data point i and the nearest centroid, i.e the centroid to the specific cluster, and b(i) is the distance to the second nearest centroid. This will give a measure of how close a certain point is to switch to another cluster than the one it is assigned. By then taking the mean of all s-values for each k-value a measure of how good the k-means clustering is obtained. In figure 7 is the mean of all s-values plotted for each k-value.

Figure 7: The mean of s from equation 6 indicates how strong the clustering become for each k-value.

According to the local maxima in figure 7, k = 8 was chosen. There is mainly two problems with choosing the k-value empirically. First, a too low k-value leads to that some aspects in the brain dynamics can be missed. Second, a too high k-value can result in that the magnitude of connectivity change sign which leading to a erroneous correlation, this is sometimes called oversplitting [4]. However, a k-value of 8 is a reasonable choice since this number is often found in the neuroscience literature, for example by Thompson and Fransson [4] and by Allen et al [27]. This also makes the results more comparable to previous work. It was also verified that all the clusters consisted of almost the same number of time-points, this to make sure that no cluster was almost empty and that the time-points were distributed over all subjects. Finally, also other k-values were tested to see the dependence of the choice of k.

(24)

3.10 State-Graphlets and Temporal-Graphlets

To correlating time-points at each cluster between each ROI-combination (edge), the Spearman correlation was used. This was used since it is robust to potential outliers that might occur due to the clustering process. For the Power 264 atlas, this resulted in one 264×264 connectivity matrix (106×106 for the HO-atlas) for each cluster with a specific connectivity value, ρ for each edge. These matrices are called state-graphlets or s-graphlets. In figure 8 A-B, an example of the step from a PCA time-series to the creation of a s-graphlets are visualized.

By representing every time-point with their respective s-graphlets in a time- series of connectivity matrices a temporal-graphlets or t-graphlets is created [4]. Thus the dynamics of the resting-state networks such as burstiness can be computed on the t-graphlets. In figure 8, an example of the construction of s-graphlets and t-graphlets is shown for a two ROI network.

(25)

Figure 8: Visualization of the construction of s-graphelts and t-graphlets for a two ROI network. A) The two ROIs and their time-series are shown, the colorbar shows the cluster affiliation obtained by the k-means clustering. B) The correlation, ρ between the two nodes for each of the clusters is then representing the edge value in the s-graphlet for these nodes. C) The time-series for the edge obtained when constructing the t-graphlet by putting the s-graphlets back to a time-series. (Image source: Bursty and persistent properties of a large- scale brain networks revealed with a point-based method for dynamic functional connectivity [4]).

3.11 Weighting

The t-graphlets now only consists of eight unique s-graphlets in a time-series which means that it only consists of eight unique connectivity values for each

(26)

edge. This is not ideal as we would like unique connectivity estimates per time- point. To get a unique value for each edge at each time-point the t-graphlets are weighted. To do this, two different weighting strategies have been developed.

Both weighting schemes aims to consider two variables, the s-graphlet and the distance from the cluster centroid. Given these considerations, the time-series shall be manipulated to consist of unique connectivity values for each edge. By considering each s-graphlet as representative of the connectivity for time-points clustered to that centroid and assume that points closer to the centroid reflect the connectivity better than points further away. This assumes that there is some linearity between the distance to a cluster’s centroid and how much the connectivity is represented by the said cluster. The two weighting schemes starting from the same equation 7

T Gi,j,t=

k

X

s=1

SGi,j,swi,j,s,t (7)

which is simply that the edge, i, j at the time, t for a t-graphlet, T G is the weighted mean of all the s-graphlets, SG were the weight, w is specific for the cluster, s and time-point, t. The difference between the two methods is of course the weight which is described below.

3.11.1 Global time-point scaling

In the first method called global time-point scaling (GTPS), the weight is derived from the distance to each centroid as:

w_i,j,s,t= d^−α_s,t

k

P

s

d^−α_s,t

(8)

where d_s,tis the euclidean distance from a time-point to a centroid of cluster, s at the time point, t. Since d is identical for every edge in a specific s-graphlet, the weight, w will be a scalar applied to every edge and α is some scaling value, in this study α = 4 was used. A higher value of α will increase the significance of nearer clusters and when α increase to infinity, the value will become the closest unweighted t-graphlet. This weighting scheme is very quick and easy to implement but it effectively scales all edges of a t-graphlet by an equal amount.

This seem unreasonable since every edge may display its peak of connectivity at a specific time-point and having every connection increase or decrease based on the time-point does not seem biologically plausible.

3.11.2 Edge time-point scaling

The second weighting scheme is called edge time-point scaling (ETPS) and is similar to the first one in the idea of constructing weights derived from s- graphlets and distance from centroids. However this method aims to derive weights for each edge, instead of each s-graphlet.

The distance from each centroid is expressed along 104 PCA dimensions, mean- ing some time-points with the same Euclidean distance from a cluster in the previous method could vary along different PCA dimensions. Some of these PCA dimensions will then have a greater importance for the edge. For each

(27)

time-point and each edge, this approach considers to which extent, the PCA dimensions that the clustering was based upon explain the variance for the nodes the edge connects. This builds upon the idea that if a time-point is identical in distance from two centroids except for one PCA dimension, but this PCA dimension explains none of the variance for two nodes in a s-graphlet, then the weights for this edge should be identical. Distance for each edge becomes based on the distance of the relevant PCA components for a nodes-pairing, giving unique scalars, per time-point, per edge.

PCA express the original data, y, in PCA dimensions, x and loadings l, which is used in this weighting scheme. The loadings represent the weight to transform the PCA dimensions back to the original spatial dimensions. Hence can the original data be obtained by the following formula

y = xl (9)

First, the eucldiean distance to a cluster’s centroid, C for each PCA dimension is calculated, instead of the shortest distance used in GTPS. This is described in the following equation

dc,s,t= |Cc,s− xc,t| (10)

Where dc,s,t is the distance in each of the 104 PCA-component, c, for each cluster, s and for every time-point, t. Cc,s is the centroid and xc,t is the data expressed in the PCA dimensions. But since this method aims to obtain a specific weight for each edge is the loadings for the nodes which is connected by the current edge summed for every state as can be seen in equation 11 below

l^e_c,i,j= lc,i+ lc,j (11)

Where l^e_c,i,j is the loading in each PCA dimension, c for the edge between the nodes, i and j. For simplicity negative loadings are set to zero. Otherwise a correlation at a s-graphlet may effect a correlation negatively, this is described in the equation below.

l_c,i,j^e =

®l^e_c,i,j if l^e_c,i,j> 0

0 if l^e_c,i,j≤ 0 (12)

To obtain the weight is the edge-loadings and the distance multiplied and transformed into percent. Finally are also the weight scaled into percent in the same way as in equation 8. This is described in equation 13.

wi,j,s,t=

Ç _d_c,s,t_l^e

c,i,j

P

c

dc,s,tl^e_c,i,j

å⁻¹

k

P

s

Ç d_c,s,tl^e_c,i,j

P

s

d_c,s,tl^e_c,i,j

å⁻¹ (13)

Where wi,j,c,t is the final weight for the edge between the nodes, i and j, for every cluster, s and for all the time-points, t. This is then used in equation 7 to obtain the weighted t-graphlet.

3.12 Thresholding

Since it is unreasonable to assume that every region of the brain is connected at every time-point is the t-graphlets thresholded and since only the number of

(28)

bursts and not the magnitude of the connectivity during a burst is of interest the t-graphlets are made binary. For the binary t-graphlets are 1 representing a present connectivity at the edge and 0 representing a absent connectivity at the edge.

In the article The mean– variance relationship reveals two possible strategies for dynamic brain connectivity analysis in fMRI [28] a negative correlation between the mean and the variance for each edge time-serie is suggested. This lead too mainly two possible types of thresholds. These are variance threshold and magnitude threshold. Since the connectivity is higher within resting state networks also the mean will be higher, this implies that the variance will be lower according to the negative scaling [28]. This lead to that the variance based threshold is more sensitive to the between resting state network connectivity. While the magnitude threshold is more sensitive to the within resting state network connectivity [28]. Considering this projects aims to calculate the integration of information across the entire brain, the magnitude threshold was not used, as it would results in far more bursts within network nodes compared to between network connectivity. This would lead the method to be biased by the size of each resting-state network, which is not what the goal is. Instead was the variance threshold used, which is described in section 3.12.1. Also a new threshold, the global variance threshold was developed, this is described in section 3.12.2. But before thresholding the data is demeaned along the time dimension and subjects dimension, so that the connectivity will be around zero for all time-series. This is done to make sure that every node will be treated independent of the magnitude of the mean, why this is important will soon be discussed.

3.12.1 Local Variance Threshold

The local variance threshold (LVT) is very straight forward. The edge is rejected (i.e set to zero) if the connectivity value is lower then two standard deviations for the current edges time-series. And if it is over two standard deviations it is set to 1. This thresholding strategy entails that every edge will have some activity that excedes the threshold. This may add noisy connections but at the same time does not bias the edge selection. In figure 9 is an example showing how the LVT thresholds the time-series.

(29)

Figure 9: Example of how the LVT thresholds two different time-series.

3.12.2 Global Variance Threshold

The global variance threshold (GVT) is more complicated. First a so called null- s-graphlets is created. The difference between the null-s-graphlet and the original s-graphlet is that one of the dimensions order of the time-points clustered in a s-graphlet was flipped. This means that the null-s-graphlets mainly consists of noise and any effects due to an autocorrelation. Next, the null-s-graphlets go through the same step as the s-graphlets to become a null-t-graphlet i.e weighting and demeaning. Every edge from the null-t-graphlets are combined to create a distribution of correlations which we can compare the t-graphlet values against to see if the t-graphlet correlations are large enough in magnitude to be considered relevant. The threshold is then set so that the edge is accepted if its connectivity is higher than 99% (p < 0.01) of the edges in the null-t-graphlet.

Since this will end up in just a single number, the demeaning is very important, without this would this threshold become a magnitude threshold. An example of the GVT can be seen in figure 10. This thresholding attempt reduces the number of noisy edges having time-points marked as relevant connectivity.

However, this will lead to the edges with the highest variance to have more time-points marked as relevant (see figure 11E).

(30)

Figure 10: Example of how the GVT thresholds two different time-series.

3.13 Bursty centrality

With bursts being increases in connectivity under a short time-period. Bursty centrality aims to, for every node, quantify how many bursts occur throughout the time-series. To do this, the number of bursts exists within a connectivity time-series needs to be determined. For doing this automatically, two parameters is needed: one define the minimum length of a burst (min_L), i.e the number of ones that a burst must include. And one, inter spike interval (ISI), which define how many consecutive time-points with zero connectivity are allowed between connectivity time-points to still be considered a burst.

In this project mainly ISI = 2 and minL= 1 was used, but also other values of these parameters were tested and compared, for example ISI = 0 and minL = 1 (figure 18). By the use of these settings could now the number of bursts for every edge be calculated, this is defined as counts. By then taking the mean of this for every node is the bursty centrality, BC obtained, this is described in equation 14

BCi= 1 n

X

j

Ci,j (14)

Where Ci,j is the number of bursts in the edge between the nodes, i, j and n is the number of nodes. To make it easier to compare bursty centrality to other graph theoretical measures, the bursty centrality is scaled between 0 and 1, this can be seen in the equation below.

BC_i^∗= BCi− min[BCi]

max[BC_i] − min[BC_i] (15)

Where BC_i^∗ is the scaled bursty centrality for each node, i.

(31)

Dynamic fMRI brain connectivity: 4 Results

4 Results

As can be seen in the method in section 3, the methods robustness is challenged for several different options, for example the choice of atlas (i.e Power 264 or Harvard-Oxford), weighting method (i.e the GTPS or ETPS) and k-value. But there is also a choice of thresholding algorithm that have to be customized according to the present study. The dependence of all these are studied and is presented in this section. The method was also tested for a new dataset to see whether there is any common patterns or common nodes that are the most bursty across these datasets. Finally was also the bursty centrality correlated with other graph theoretical measures such as betweenness centrality and temporal centrality. If nothing else is stated, the different parameters were computed with Oulu A dataset, Power 264 atlas parcellation, ETPS and with k = 8.

4.1 Dependence of thresholding algorithm

To replicate the findings in The mean– variance relationship reveals two possible strategies for dynamic brain connectivity analysis in fMRI [28], which was discussed in section 3.12 was the mean plotted versus the variance (figure 11A).

This plot was the foundation for the two different thresholding algorithms. How- ever, in this case, there is no linear correlation between the mean and the variance, but a similar dissociation between edges with high mean and high variance exists. Still it is important to understand how the different thresholding algorithms thresholds the data. For both the LVT and GVT algorithms is the number of counts (i.e number of bursts in every edge) almost the same regardless of the mean (figure 11B for LVT and figure 11D for GVT). However, in the GVT algorithm a clear relationship can be seen between the counts and the variance (figure 11E) which does not exist with the mean (figure 11D). Since GVT uses the global variance as a threshold for every edge is this a quite obvious behavior.

(32)

Figure 11: A) Mean versus variance of the edges time-series. B) Shows that there is no clear dependence between the number of counts (i.e number of bursts in every edge) for the LVT and the mean for the edge time-series. C) Shows that there is no clear dependence between the number of counts for the LVT and the variance for the edge time-series. D) Shows that there is no clear dependence between the number of counts for the GVT and the mean for the edge time- series. E) Shows a clear positive dependence between the number of counts and variance of the edge time-series when using the GVT.

We have established that the different thresholding algorithms will allow different edges to pass the threshold and it is of interest to see which nodes in the brain that is the most bursty. The bursty centrality was calculated for each node (figure 12A for the LVT and figure 12B for the GTV). The nodes with the highest bursty centrality were classed as candidates to be dynamic hubs.

A threshold at 15% (40 nodes) was chosen, which is the yellow part in figure 12 A and B and the colorbar along the x-axis show the network corresponding to each node. These 40 nodes corresponds to approximately 23% of the total number of bursts for the LVT and 32% for the GVT.

Using BrainNet viewer were these 40 nodes plotted for both the LVT (figure 12C)) and for the GVT (figure 12D). The size of the ROIs corresponds to the number of bursts and the color corresponds to the network affiliation. These figures shows that for the LVT, the most bursty nodes are located in the frontal parts of the brain, the most bursty node is left thalamus in the subcortical network. In table 1 in appendix A can a complete list of the 40 most bursty nodes be found. Worth mentioning is that the DMN is by far the most bursty network where 14 of the 40 most bursty nodes are located. For the GVT on the other hand is there mainly the dorsal regions that are bursty, the most bursty one is left precuneus in the somato-motor network. For this threshold the visual and the DMN that are dominating with 13 respective 11 nodes of the 40 most bursty. The complete list of this can be found in table 2 in appendix A.

(33)

Figure 12: A-B) Bursty centrality for the (A) LVT and (B) GVT. The colorbar shows the network affiliation and the yellow part shows the 40 nodes (15%) with highest bursty centrality. C-D) The 40 nodes with highest bursty centrality plotted using BrainNet viewer, for (C) LVT and (D) GVT. The size of the nodes indicates the bursty centrality and the color shows the network affiliation

(34)

4.2 Comparison between atlases

To investigate how the result will depend on the choice of atlas, the plots corresponding to figure 11 and 12 A and B were also plotted for the HO-atlas described in section 3.7, these can be seen in figure 13. Since the network clas- sification in the HO-atlas only are based on mutual information and not further investigated as in the case for the Power 264 atlas is there no colorbar in figure 13 F and G. Instead, the nodes are listed in table 3 for the LVT in figure 13F and table 4 for the GVT in figure 13G, these tables can be found in appendix B. By comparing figure 13 to figure 11 can it be observed that the mean versus variance plot is very similar for the both atlases. It can also be seen that the thresholds treats the data in the same way as previously. Also the distribution of bursts following the same shape in both cases, for the HO-atlas corresponds the 15% (i.e 16 nodes) most bursty nodes for approximately 17% of the total number of bursts for the LVT, where left heschl’s gyrus have the highest bursty centrality. For the GVT have frontal medial cortex the highest bursty centrality and the 16 most bursty nodes corresponds for approximately 31% of the total number of bursts.

(35)

Figure 13: A) The mean versus variance of the edges time-series plotted for the HO-atlas, this shows the same behavior as for the Power 264 atlas. B) Shows that there is no clear dependence between the number of counts (i.e number of bursts in every edge) for the LVT and the mean for the edge time-series. C) Shows that there is no clear dependence between the number of counts for the LVT and the variance for the edge time-series. D) Shows that there is no clear dependence between the number of counts for the GVT and the mean for the edge time-series. E) Shows a clear positive dependence between the number of counts and variance of the edge time-series when using the GVT for the HO- atlas. F-G) Bursty centrality for the (F) LVT and (G) GVT, the yellow part shows the 16 (15%) nodes with highest bursty centrality. The nodes order can be seen in table 3 (LVT) respectively 4 (GVT) in appendix B

4.3 Dependence of k-value

Further, also the dependence of the k-value was investigated. To do this, the bursty centrality was calculated for all nodes with k-values ranging from 5 to 14. The correlation between these was then plotted into a correlation matrix which can be seen in figure 14A for LVT and 14B for GVT. For the LVT is the correlation coefficient, ρ ranging from 0.1579 between k = 7 and k = 11 up to 0.6858 between k = 12 and k = 13. For the GVT in figure 14B is there considerably higher correlation. In this case is the correlation coefficient ranging from 0.7302 between k = 5 and k = 8 up to 0.9464 between k = 10 and k = 11.

Noteworthy is that the correlation coefficient seems to be considerably lower for the lowest number of k, especially for k = 5 to k = 8 for the LVT and k = 5 for the GVT. So for the both thresholds is there a correlation that is significantly different from zero, however the correlation is considerably much stronger when using the GVT.