# Resting state FSLnets practical

In this practical you will perform network modelling using fMRI timeseries data (i.e. the output from stage 1 of dual regression). The toolbox you will use for this is FSLNets, which currently requires MATLAB (or Octave - a free alternative).

## Contents:

Networks estimation
Estimating network matrices from dual regression outputs
Group-average netmat summaries
Calculating group-average netmats
Cross-subject comparison with netmats
Comparing individual edge strengths between subject groups
Multivariate cross-subject analysis (optional)
Multivariate comparison of whole netmats across subject groups

Estimating network matrices and performing statistical analyses on these requires the following steps:

1. Extract subject-specific timecourses relating to a given set of spatial maps of nodes. For example, using dual regression stage 1 outputs
2. Look for group-level noise components that you want to remove
3. Calculate correlations (full or partial) between all pairs of timeseries to build the network matrix (netmat) for each subject
4. Look at the group-average network matrix and how nodes cluster together into a hierarchy
5. Perform group-level statistical analysis

In the next sections we will go through an example of how to run these steps. For this practical we re-ran for you the group ICA and the dual regression using 100 components, with the following command (please do NOT run this again!):

```melodic -i input_files.txt -o groupICA100 \
--tr=0.72 --nobet -a concat \
--bgimage=\$FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz \
--report --Oall -d 100
```

## Networks estimation

To start this FSLNets network modelling practical, open a terminal and cd to the working directory of this practical, and then start octave:

```cd ~/fsl_course_data/rest/Nets
octave
```

Add the folder containing the scripts we are going to use to the path, and setup a few variables:

```addpath ~/fsl_course_data/rest/FSLNets
```

Now we are ready to start setting up some important parts for our network modeling analysis. To load in all subjects' timeseries data files from the dual-regression output directory, run:

```ts = nets_load('groupICA100.dr', 0.72, 0);
```

Type fieldnames(ts) and hit enter to have a look at the contents of the structure that was created using this command. ts.ts contains the stage 1 dual regression timeseries from all subjects. Try to figure out what information the other variables in the structure ts contain.

digraph G { rankdir=LR; node [shape=box]; 1 [label="dual regression output folder ('groupICA100.dr')"]; 2 [label="repetition time (TR) (0.72)"]; 3 [label="normalisation of timeseries (0) here none"]; 4 [label="nets_load"]; 5 [label="output structure with subfields (ts)"]; 1 -> 4; 2 -> 4; 3 -> 4; 4 -> 5; }

We will now have a quick look at the temporal spectra of the RSNs, as a quick check that these look reasonable. Running the following command will show you a figure in which the left part of the plot shows one spectrum per group-ICA component, each averaged across all subjects. The right part shows the same thing, but with all spectra overlapping.

```ts_spectra = nets_spectra(ts);
```

Why do the figures show that all nodes have high power in the low-frequency range?
Incorrect! While the TR does relate to the maximum frequency that can be measured, it does not determine where the power of the signal is.
Incorrect! These data have only been highpass filtered to remove drift.
Correct! Neural responses are measured indirectly in BOLD, and fluctuations are expected to vary slowly as a result of the properties of the HRF.

There is now the option to remove components' (nodes') timeseries that correspond to artefacts rather than plausible nodes. Similar to the last practical, you need to decide this by looking at the spatial maps, timecourses, and frequency spectra. To save time, we have listed the good components for you. Run the following commands to list the good components and apply the cleanup:

```ts.DD = [1:3,5,6:9,11:13,17:23,25:38,40,42,43,47:50,52,53,55:59,61,...
62,64:66,70:74,77,80,81,86,87,93,97];
ts = nets_tsclean(ts,1);
```
digraph G { rankdir=LR; node [shape=box]; 1 [label="input structure (ts) loaded previously"]; 2 [label="aggressive or soft cleanup (1) i.e. aggressive"]; 3 [label="nets_tsclean"]; 4 [label="output structure (ts) overwriting original ts"]; 1 -> 3; 2 -> 3; 3 -> 4; }

What is 'aggressive' cleanup
Correct!
Incorrect! Simply removing timeseries is called soft cleanup.

Now you are ready to compute a network matrix for each subject, which is in general a matrix of connection strengths. We will compute two different versions of these netmats for each subject: a simple full correlation ('corr'), and a partial correlation that has been regularised in order to potentially improve the mathematical robustness of the estimation ('ridgep'). The partial correlation matrix should do a better job of only estimating the direct network connections than the full correlation does.

```Fnetmats = nets_netmats(ts,1,'corr');
Pnetmats = nets_netmats(ts,1,'ridgep',0.1);
```
digraph G { rankdir=LR; node [shape=box]; 1 [label="inputs structure (ts)"]; 2 [label="apply Fisher's r-to-z transformation? (1) i.e. yes"]; 3 [label="method for netmat estimation ('ridgep') i.e. regularised partial correlation"]; 4 [label="regularisation parameter (0.1)"]; 5 [label="nets_netmats"]; 6 [label="output variable (Pnetmats)"]; 1 -> 5; 2 -> 5; 3 -> 5; 4 -> 5; 5 -> 6; }

The full and partial netmats are now calculated for all subjects. Now run this command to show you the size of the Fnetmats variable:

```size(Fnetmats)
```

Can you figure out how the netmats of all subject are saved in this? Answer. In the next section, you are going to compute the group average and take a look at the average network matrix.

## Group-average netmat summaries

Now that we have computed the full and partial network matrices for each individual subject. The next step is to perform a simple group-level analysis to look at the mean network matrix across all subjects. This can be done using the command below, which saves out both the simple average (Mnet) and the results of a simple one-group t-test across subjects as Z values (Znet).

```[Znet_F,Mnet_F]=nets_groupmean(Fnetmats,0);
[Znet_P,Mnet_P]=nets_groupmean(Pnetmats,1);
```

The second input to this command decides whether or not to display a summery figure. This is set to 1 for the second command run on the partial netmats, so the summary figure for this will show up. On the left, this figure shows the results from the group t-test, and on the right is a consistency scatter plot showing how similar the results from each subject are to the group (i.e. the more this looks like a diagonal line, the more consistent the relevant netmat is across subjects).

The value of `Mnet_P` in row 3, column 27 is ~6.6 (you can check this by typing `Mnet_P(3,27)` into Octave and hitting enter). What does this number represent?
Incorrect! Remember that we are visualising the results for the partial netmats!
Correct! The group-averaged strength of the partial correlation between the pair of nodes is shown here (which is also called the 'edge' between the two nodes). The number is higher than 1, because these are averaged z-transformed correlation coefficients.
Incorrect! We have not yet done a group comparison, this is just the average partial correlation strength between the pair of nodes (which is called an 'edge').

The next thing we can look at is how nodes cluster together to form larger resting state networks. For this we run a clustering method that groups nodes together based on their timeseries. To view this network hierarchy, run:

```nets_hierarchy(Znet_F,Znet_P,ts.DD,'groupICA100.sum');
```

digraph G { rankdir=LR node [shape=box]; 1 [label="group-averaged netmat to drive clustering (Znet_F) shown below the diagonal"]; 2 [label="group-average netmat to show (doesn't drive clustering) (Znet_P) shown above the diagonal"]; 3 [label="list of good components (ts.DD) entered earlier"]; 4 [label="directory containing png summary figures ('groupICA100.sum') to create run: slices_summary groupICA100/melodic_IC 4 \$FSLDIR/data/standard/MNI152_T1_2mm groupICA100.sum -1"]; 5 [label="nets_hierarchy = the FSLnets command we are running"]; 1 -> 5; 2 -> 5; 3 -> 5; 4 -> 5; }

You can see, for example, that the nodes group together in the dark blue tree on the far left are part of a large-scale resting state network called the default mode network that you may have heard about.

## Cross-subject comparison with netmats

We are now able to test whether the netmats differ significantly between healthy controls and patients with a tumor using a two-sample t-test. This is a 'univariate' test, as we will test each network matrix edge separately for a group-difference, and then we will estimate p-values for these tests, correcting for the multiple comparisons across all edges. By analogy to high-level task-fMRI analyses: you can think of each subject's netmat as being an NxN image of voxels, and the univariate testing as modelling each voxel (in isolation from each other) across subjects. To perform the comparison please run (which calls `randomise` from within Octave, with 5000 permutations for each contrast):

```[p_uncorr,p_corr]=nets_glm(Pnetmats,'design/unpaired_ttest_1con.mat','design/unpaired_ttest_1con.con',1);
```
digraph G { rankdir=LR; node [shape=box]; 1 [label="input to t-test (Pnetmats)"]; 2 [label="EVs of group-level design ('design/paired_ttest.mat') use GLM gui to create this"]; 3 [label="contrasts of group-level design ('design/paired_ttest.con') use GLM gui to create this"]; 4 [label="create output figure? (1) i.e. yes"]; 5 [label="nets_glm"]; 6 [label="variable containing uncorrected p-values (p_uncorr)"]; 7 [label="variable containing p-values corrected for multiple comparisons (p_corr)"] 1 -> 5; 2 -> 5; 3 -> 5; 4 -> 5; 5 -> 6; 5 -> 7; }

Once randomise has finished, you will see a figure showing "netmats" containing corrected p-values. The results above the diagonal show edges where the two-group t-test is significant, at corrected-p<0.05.

We will now run a command that shows which nodes were linked to the significant result:

```nets_edgepics(ts,'groupICA100.sum',Znet_P,reshape(p_corr,ts.Nnodes,ts.Nnodes),1);
```

In addition, we also want to show how the partial correlation differs between the patients and the controls and these two significant edges. To do this, run:

```nets_boxplots(ts,Pnetmats,57,33,6);
```

The boxplots summarize the distributions of the correlation values (connections strengths) in the two groups - A being healthy controls and B being tumour patients - for this one particular node-pair (57,33).

Using this plot, how would you interpret what caused the significant group difference for the edge between node 57 and 33?
Correct! For the patients the results are not significantly different from zero.
Incorrect! The group-averaged strength of the partial correlation was positive in healthy controls, and close to zero in patients.
Incorrect! The group-averaged strength of the partial correlation was positive in healthy controls, and close to zero in patients.

Type `exit` in the terminal to exit octave if you are finishing the practical here and not doing the optional section below.

## Multivariate cross-subject analysis (optional)

Finally, we will now carry out multivariate cross-subject analysis - meaning that, instead of considering each netmat edge in isolation (like we did above) we will consider the whole netmat in one analysis.

For example, we can attempt to classify subjects into patients or controls using machine learning methods, such as support vector machines (SVM) or linear discriminant analysis (LDA). Such methods look at the overall pattern of values in the netmat, and try to learn how the overall pattern changes between the two groups.

The following command feeds the regularised partial correlation netmats from both groups into LDA. It uses a method known as leave-one-out cross-validation to train and test a classifier, and reports in what percentage of tests it was successful at discriminating between patients and controls:

```nets_lda(Pnetmats,6,2)
```

You can also try using other classifiers, for example by changing the last entry from 2 to 1 (type help nets_lda for further details).

One "downside" of such multivariate testing is that you can no longer make strong statistical claims about individual edges in the network - the whole pattern of edges has been used, so we don't know which individual edges are significantly different in the two groups.

Type `exit` in the terminal to exit octave.

The End.