Model

For details on the model, please see SIEM model.

Real Data Experiments

In this notebook, we determine whether there exists a difference in connectivity \(p_{homo}\) homotopically (same region, opposite hemisphere) vs \(p_{hetero}\) heterotopically (different region) connectivity within a particular modality. We consider \(\delta_{x} = p_{homo} - p_{hetero}\) to be the difference in connectivity for a graph from a graph or collection of graphs from a particular modality \(x\).

Test

Our test for this notebook is as follows:

\[\begin{align*} H_0: \delta_d \leq \delta_f \\ H_A: \delta_d > \delta_f \end{align*}\]

in words, whether the connectivity homotopically exceeds the connectivity heterotopically. We will do this in 2 ways:

  • An \(n\) sample test to determine whether the difference in connectivity homotopically vs heterotopically in the functional connectomes exceeds that of the diffusion connectomes at the population level.
  • A \(2\)-sample test to determine whether the difference in connectivity homotopically vs heterotopically in the functional connectomes exceeds that of thet diffusion connectomes at the individual graph level (1 diffusion and 1 functional graph from the same subject).

We will perform this experiment for both dMRI and fMRI-derived connectomes.

Raw Data

For the data, we compute the weighted mean functional (rank of each edge) and diffusion (number of fibers). For the functional connectome, we threshold such that the largest 50% of edges are set to connected, and the smallest 50% set to disconnected. For the diffusion (which are natively sparse) we just threshold edges that are present to connected, and edges that are not present to disconnected (threshold about 0).

The data below can be downloaded and moved to appropriate folders as follows (note that the below section requires sudo access):

sudo mkdir /data/
sudo chmod -R 777 /data

cd /data
wget http://openconnecto.me/mrdata/share/derivatives/dwi_edgelists.tar.gz
wget http://openconnecto.me/mrdata/share/derivatives/fmri_edgelists.tar.gz
wget http://openconnecto.me/mrdata/share/connectome_stats/connectome_stats.zip

mkdir -p /data/connectome_stats /data/all_mr /data/all_mr/dwi/edgelists /data/all_mr/fmri/ranked/edgelists
mv dwi_edgelists.tar.gz /data/dwi/edgelists
cd /data/dwi/edgelists
tar -xvzf dwi_edgelists.tar.gz
mv /data/fmri_edgelists.tar.gz /data/fmri/ranked/edgelists
cd /data/fmri/ranked/edgelists
tar -xvzf fmri_edgelists.tar.gz
mv /data/connectome_stats.zip /data/connectome_stats.zip
cd /data/connectome_stats
unzip connectome_stats.zip
basepath = '/data/connectome_stats/'
fmri_gr = read_graph(file.path(basepath, 'fmrimean_1709.edgelist'), format="ncol")
vset <- V(fmri_gr)
ordered_v <- order(vset)
fmri_gr = read_graph(file.path(basepath, 'fmrimean_1709.edgelist'), format="ncol", predef=ordered_v)
fmri_mean = get.adjacency(fmri_gr, type="both", sparse=FALSE, attr='weight')
dwi_gr = read_graph(file.path(basepath, 'dwimean_2861.edgelist'), format="ncol", predef=ordered_v)
dwi_mean = get.adjacency(dwi_gr, type="both", sparse=FALSE, attr='weight')

fmri_thresh = thresh_matrix(fmri_mean)
dwi_thresh = thresh_matrix(dwi_mean, thresh=0)

gs.plot.plot_matrix(fmri_thresh, title = "Mean Thresholded Functional Connectome", legend.name = "connection", ffactor = TRUE) +
  theme_bw()

gs.plot.plot_matrix(dwi_thresh, title = "Mean Thresholded DWI Connectome", legend.name = "connection", ffactor = TRUE) +
  theme_bw()

Blocked Data

here, we will compute the probability of an edge existing in each of 4 quadrants (2 homotopic quadrants; 2 heterotopic quadrants):

group1 = 1:35
group2 = 36:70
groups = list(group1, group2)
fmri_block = block_data(fmri_thresh, groups)
dwi_block = block_data(dwi_thresh, groups)

fmriu.plot.plot_graph(fmri_block, title = "Blocked Functional Connectome", xlabel = "Hemisphere",
                      ylabel="Hemisphere", include_diag = TRUE, legend.name = "p")

fmriu.plot.plot_graph(dwi_block, title = "Blocked DWI Connectome", xlabel = "Hemisphere",
                      ylabel="Hemisphere", include_diag = TRUE, legend.name = "p")

Diffusion

nroi <- 70
dwi.dsets = c('BNU1', 'BNU3', 'HNU1', 'KKI2009', 'NKI1', 'NKIENH', 'MRN1313', 'Templeton114', 'Templeton255', 'SWU4')
dwi.atlas = 'desikan'
dwi.basepath = '/data/all_mr/dwi/edgelists'

graphobj = fmriu.io.collection.open_graphs(basepath = dwi.basepath, atlases = dwi.atlas, datasets = dwi.dsets,
                                           gname = 'graphs', fmt='edgelist', rtype = 'array')
## [1] "opening graphs for BNU1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for BNU3 dataset and desikan parcellation atlas..."
## [1] "opening graphs for HNU1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for KKI2009 dataset and desikan parcellation atlas..."
## [1] "opening graphs for NKI1 dataset and desikan parcellation atlas..."
## [1] "opening graphs for NKIENH dataset and desikan parcellation atlas..."
## [1] "opening graphs for MRN1313 dataset and desikan parcellation atlas..."
## [1] "opening graphs for Templeton114 dataset and desikan parcellation atlas..."
## [1] "opening graphs for Templeton255 dataset and desikan parcellation atlas..."
## [1] "opening graphs for SWU4 dataset and desikan parcellation atlas..."
dwi.graphs = graphobj$graphs
dwi.datasets = graphobj$dataset
dwi.subjects = graphobj$subjects
ne = 1225
nroi <- 70
group1 <- c()  # edges in same hemisphere
group2 <- c()  # edges across hemispheres
for (i in 1:nroi) {
  for (j in 1:nroi) {
    idx <- (i - 1)*nroi + j
    if (abs(j - i) == 35) {  # across hemispheric edges
      group1 <- c(group1, idx)
    } else if (i != j) {  # ignore diagonal
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
dwi.models <- sapply(1:dim(dwi.graphs)[1], function(i) {
                  gs.siem.fit(thresh_matrix(dwi.graphs[i,,], 0), Es, alt='greater')
                }, simplify = FALSE)

dwi.homo.phat <- sapply(dwi.models, function(model) model$pr[1])
dwi.hetero.phat <- sapply(dwi.models, function(model)model$pr[2])
dwi.delta.phat <- sapply(dwi.models, function(model) model$dpr[1,2])
dwi.delta.var <- sapply(dwi.models, function(model)model$dvar[1,2])
dwi.phat.mu = mean(dwi.delta.phat)
dwi.phat.var = model.var(mean(dwi.homo.phat), ne) + model.var(mean(dwi.hetero.phat), ne)

Functional

nroi <- 70
fmri.dsets = c('BNU1', 'BNU2', 'BNU3', 'HNU1', 'IBATRT', 'IPCAS1', 'IPCAS2', 'IPCAS5', 'IPCAS6', 'IPCAS8', 'MRN1', 'NYU1', 'SWU1', 'SWU2', 'SWU3', 'SWU4', 'UWM', 'XHCUMS')
fmri.atlas = 'desikan-2mm'
fmri.basepath = '/data/all_mr/fmri/ranked/edgelists/'

graphobj = fmriu.io.collection.open_graphs(basepath = fmri.basepath, atlases = fmri.atlas, datasets=fmri.dsets, fmt='edgelist', rtype = 'array')
## [1] "opening graphs for BNU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for BNU2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for BNU3 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for HNU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IBATRT dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS5 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS6 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for IPCAS8 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for MRN1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for NYU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU1 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU2 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU3 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for SWU4 dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for UWM dataset and desikan-2mm parcellation atlas..."
## [1] "opening graphs for XHCUMS dataset and desikan-2mm parcellation atlas..."
fmri.graphs = graphobj$graphs
fmri.datasets = graphobj$dataset
fmri.subjects = graphobj$subjects
ne = 1225
nroi <- 70
group1 <- c()  # edges in same hemisphere
group2 <- c()  # edges across hemispheres
for (i in 1:nroi) {
  for (j in 1:nroi) {
    idx <- (i - 1)*nroi + j
    if (abs(j - i) == 35) {  # across hemispheric edges
      group1 <- c(group1, idx)
    } else if (i != j) {  # ignore diagonal
      group2 <- c(group2, idx)
    }
  }
}
Es <- list(group1, group2)
fmri.models <- sapply(1:dim(fmri.graphs)[1], function(i) {
                  gs.siem.fit(thresh_matrix(fmri.graphs[i,,], 0.5), Es, alt='greater')
                }, simplify = FALSE)

fmri.homo.phat <- sapply(fmri.models, function(model) model$pr[1])
fmri.hetero.phat <- sapply(fmri.models, function(model)model$pr[2])
fmri.delta.phat <- sapply(fmri.models, function(model) model$dpr[1,2])
fmri.delta.var <- sapply(fmri.models, function(model)model$dvar[1,2])
fmri.phat.mu = mean(fmri.delta.phat)
fmri.phat.var = model.var(mean(fmri.homo.phat), ne) + model.var(mean(fmri.hetero.phat), ne)

\(n\) sample test

Here, we take each functional and diffusion connectome individually, and compute the parameters of our block model for each connectome. The question we seek to first answer is, given a large number of observations of \(\hat{\delta}\), can we detect a difference in homotopic vs heterotopic connectivity in the diffusion connectomes compared to the functional connectomes?

We might want to visualize the distribution of \(\delta = \hat{p}_{homo} - \hat{p}_{hetero}\) under the analytical model and compare to our empirical model for functional and diffusion:

ne = 1225
# density estimates of the two populations of delta
dwi.distr.emp.mod = density(as.numeric(dwi.delta.phat))
fmri.distr.emp.mod = density(as.numeric(fmri.delta.phat))

# variances sum for analytical computation
dwi.distr.ana = dnorm(dwi.distr.emp.mod$x, mean=dwi.phat.mu, sd=sqrt(dwi.phat.var))
fmri.distr.ana = dnorm(fmri.distr.emp.mod$x, mean=fmri.phat.mu, sd=sqrt(fmri.phat.var))

n_diff = length(dwi.distr.emp.mod$x)
dwi.dat = data.frame(x = c(dwi.distr.emp.mod$x, dwi.distr.emp.mod$x), y = c(dwi.distr.emp.mod$y, dwi.distr.ana),
                      distribution=c(rep("empirical", n_diff), rep("analytical", n_diff)))
dwi.dat$distribution = factor(dwi.dat$distribution)

ggplot(dat=dwi.dat, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta_D$')) +
  ggtitle(TeX('Distribution of $\\delta_D = \\hat{p}_{homo} - \\hat{p}_{hetero}$, DWI')) +
  theme(panel.background = element_rect(fill = '#ffffff'))

n_diff = length(dwi.distr.emp.mod$x)
fmri.dat = data.frame(x = c(fmri.distr.emp.mod$x, fmri.distr.emp.mod$x), y = c(fmri.distr.emp.mod$y, fmri.distr.ana),
                      distribution=c(rep("empirical", n_diff), rep("analytical", n_diff)))
fmri.dat$distribution = factor(fmri.dat$distribution)

ggplot(dat=fmri.dat, aes(x=x, y=y, ymax=y, fill=distribution, color=distribution, group=distribution)) +
  geom_ribbon(ymin=0, alpha=0.5) +
  ylab('Density') +
  xlab(TeX('$\\delta_F$')) +
  ggtitle(TeX('Distribution of $\\delta_F = \\hat{p}_{homo} - \\hat{p}_{hetero}$, fMRI')) +
  theme(panel.background = element_rect(fill = '#ffffff'))