SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 Package Loading and Data Manipulation

library(TCIU)
library(DT)
library(R.matlab)
library(AnalyzeFMRI)
library(ggfortify)
mat1 = readMat("subj1_run1_complex_all.mat")
bigim1 = mat1$bigim[,64:1,,]
bigim1_mod = Mod(bigim1)
smoothmod = GaussSmoothArray(bigim1_mod,sigma = diag(3,3))
# dim(bigim1) = 64 64 40
# bigim1 contains the complex image space
# dimensions are 64x*64y*40z*160t, corresponding to x,y,z,time
load("mask.rda") # load the 3D nifti data of the mask
load("mask_label.rda") # load the 3D nifti data of the mask with labels
load("mask_dict.rda") # load the label index and label name
label_index = mask_dict$index
label_name = as.character(mask_dict$name)
label_mask = mask_label
load("hemody.rda") # load the hemodynamic contour of the brain

2 Time-series graphs

2.1 Interactive time-series visualization

The function fmri_time_series is used to create four interactive time series graphs for the real, imaginary, magnitude, and phase parts for the fMRI spacetime data. We can also add a reference plotly object to the plot. This function is based on the GTSplot function from package TSplotly.

2.1.1 Example fMRI(x=4, y=42, z=33, t)

reference_plot = sample[[4]]
fmri_time_series(bigim1, c(44,42,33), is.4d = TRUE, ref = reference_plot)

3 Kime-series/Kime-surfaces (spacekime analytics protocol)

This Rmd notebook uses a time-series simulation to illustrate how to transform the fMRI data at a fixed voxel location into a kime-surface (kime-series).

Notes:

  • Validate all steps in this transformation protocol and finalize this 3D visualization.

  • Try it with real fMRI data at brain voxel locations associated with the finger-tapping task.

3.1 Pseudo-code

  • Randomly generate 8 \(phi=\phi\) kime-phases for each of the 10 time radii. This yields an \(8\times 10\) array (phi_8_vec) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval \([-\pi:\pi)\), (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc.
  • Optionally, order all kime-phases (rows) from small to large for each column.
  • Fix the \(\nu=(x,y,z)\) voxel location and extract the fMRI time-course \(fMRI_{\nu}(t), \forall 1\leq t\leq 160\).
  • For binary stimuli (e.g., activation (ON) and rest (OFF) event-related design), we can split the 160-long fMRI series into 80 ON (Activation) and 80 OFF (rest) states, or sub-series.
  • Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity).
  • Convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.
  • Check for visual (graphical) and numeric differences between the fMRI intensities during the ON vs. OFF states
  • Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have a large range.
  • Generate the plot_ly text labels that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates \((-10\leq x\leq 10,-10\leq y\leq 10)\), but are computed using the polar coordinates \((1\leq t\leq 10,-\pi\leq \phi<\pi)\) and the polar-to-Cartesian transformation. The labels are \(21\times21\) arrays including the triple \((t, \phi, fMRIvalue)\). Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli).
  • Generate the \(21\times21\) kime-domain Cartesian grid by polar-transforming the polar coordinates \((t,\phi)\) into Cartesian counterparts \((x,y)\).
  • Interpolate the fMRI intensities (natively anchored at \((t,\phi)\)) to \(fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10\).
  • Use plot_ly to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain.

3.2 Function main steps illustration: Kime-surfaces / kime-series construction

3.2.1 Generate the kime-phases

Randomly generate 8 \(phi=\phi\) kime-phases for each of the 10 time radii. This yields an \(8\times 10\) array (phi_8_vec) of kime phase directions. These directions can be obtained by different strategies, e.g., (1) uniform or Laplace distribution sampling over the interval \([-\pi:\pi)\), (2) randomly sampling with/without replacement from known kime-phases obtained from similar processes, etc.

Optionally, order all kime-phases (rows) from small to large for each column.

# plot Laplacian
ggfortify::ggdistribution(extraDistr::dlaplace, seq(-pi, pi, 0.01), m=0, s=0.5)

# randomly generate 8 phi kime-phases for each of the 10 time radii
phi_8_vec <- matrix(NA, ncol=10, nrow = 8)
for (t in 1:10) { 
  # for a given t, generate 8 new phases
  set.seed(t);
  phi_8_vec[ ,t] <-
    extraDistr::rlaplace(8,mu=0,sigma=0.5)
  # rank-order the phases for consistency
  # within the same foliation leaf
  phi_8_vec[ ,t] <- sort(phi_8_vec[ ,t])
  # force phases in [-pi: pi)
  for (i in 1:8) {
    if (phi_8_vec[i,t] < -pi) 
      phi_8_vec[i,t] <- -pi
    if (phi_8_vec[i,t] >= pi) 
      phi_8_vec[i,t] <- pi
  }
}
# order all kime-phases (rows) from small to large for each column
# phi_8_vec_ordered <- apply(phi_8_vec, 2, sort)

3.2.3 Preprocess the matrices

Spatially smooth (blur) the matrices (in 2D) to reduce noise make them more representative of the process. May also need to temper the magnitude of the raw fMRI intensities, which can have a large range.

# Convert the long DF representing fMRI_ON and fMRI_OFF from polar coordinates to Cartesian coordinates
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.3-0
## Loading required package: spatstat.core
## Loading required package: nlme
## Loading required package: rpart
## spatstat.core 2.3-2
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-0
## 
## spatstat 2.2-0       (nickname: 'That's not important right now') 
## For an introduction to spatstat, type 'beginner'
matrix_ON <- matrix(0, nrow = 21, ncol = 21) 
matrix_OFF <- matrix(0, nrow = 21, ncol = 21) 
for (t in 1:10) {
  for (p in 1:8) {
    x = 11+t*cos(phi_8_vec[p,t])
    y = 11+t*sin(phi_8_vec[p,t])
    matrix_ON[x,y]  <- fMRI_ON[(p-1)*10 +t]
    matrix_OFF[x,y] <- fMRI_OFF[(p-1)*10 +t]
  }
}
# smooth/blur the matrices
matrix_ON_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_ON), sigma=0.5))
matrix_OFF_smooth <- (1/10000)*as.matrix(blur(as.im(matrix_OFF), sigma=0.5))

3.2.4 Generate plotly labels

Generate the plot_ly text labels that will be shown on mouse hover (pop-up dialogues) over each kime-surface/kime-series. These text-labels are stored in Cartesian coordinates \((-10\leq x\leq 10,-10\leq y\leq 10)\), but are computed using the polar coordinates \((1\leq t\leq 10,-\pi\leq \phi<\pi)\) and the polar-to-Cartesian transformation. The labels are \(21\times21\) arrays including the triple \((t, \phi, fMRIvalue)\). Separate text-labels are generated for each kime-surface (ONN vs. OFF stimuli).

# fix the plot_ly Text Labels
x <- vector()
y <- vector()
i <- 1
for (t in 1:10) {
  for (p in 1:8) {
    x[i] = 11+t*cos(phi_8_vec[p,t])
    y[i] = 11+t*sin(phi_8_vec[p,t])
    i <- i+1
  }
}

3.2.5 Cartesian representation

Generate the \(21\times21\) kime-domain Cartesian grid by polar-transforming the polar coordinates \((t,\phi)\) into Cartesian counterparts \((x,y)\).

hoverText <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_ON_smooth))) # tail(mytext)
custom_txt <- matrix(NA, nrow=21, ncol=21)
hoverTextOFF <- cbind(x=1:21, y=1:21, height=as.vector(t(matrix_OFF_smooth))) # tail(mytext)
custom_txtOFF <- matrix(NA, nrow=21, ncol=21)

for (x in 1:21) {
   for (y in 1:21) {
     t = sqrt((x-11)^2 + (y-11)^2)
     p = atan2(y-11, x-11)
     custom_txt[x,y] <- paste(' fMRI: ', round(hoverText[(x-1)*21+y, 3], 3),
                    '\n time: ', round(t, 0),
                    '\n phi: ', round(p, 2)) 
     custom_txtOFF[x,y] <- paste(' fMRI: ', round(hoverTextOFF[(x-1)*21+y, 3], 3),
                    '\n time: ', round(t, 0),
                    '\n phi: ', round(p, 2)) 
   }
}

3.2.6 Cartesian space interpolation

Interpolate the fMRI intensities (natively anchored at \((t,\phi)\)) to \(fMRI(x,y), \forall -11\leq x,y \leq 11, x^2+y^2\leq 10\).

xx2 <- 11 + c(-10:10) %o% cos(seq(-pi, pi, 2*pi/20))
yy2 <- 11 + c(-10:10) %o% sin(seq(-pi, pi, 2*pi/20))
#zz2 <- as.vector(matrix_ON_smooth) %o% rep(1, 21*21)
zz2 <- matrix_ON_smooth
ww2 <- matrix_OFF_smooth
dd2 <- matrix_ON_smooth-matrix_OFF_smooth

#plot 2D into 3D and make the text of the diameter (time), height (r), and phase (phi)
f <- list(family = "Courier New, monospace", size = 18, color = "black")
x <- list(title = "k1", titlefont = f)
y <- list(title = "k2", titlefont = f)
z <- list(title = "fMRI Kime-series", titlefont = f)

3.3 Function result

The following parts are the main output of our functions

3.3.1 Generate a long data-frame

Construct a data-frame with 160 rows and 4 columns; time (1:10), phases (8), states (2), and fMRI_value (Complex or Real intensity). Then, convert the long DF representing fMRI_ON and fMRI_OFF from their native (old) polar coordinates to the (new) Cartesian coordinates, using polar transformations.

datatable(fmri_kimesurface(bigim1_mod,c(44,42,33))[[1]])

3.3.2 Display kime-surfaces

Use plot_ly to display in 3D the kime-series as 2D manifolds (kime-surfaces) over the Cartesian domain.

fmri_kimesurface(bigim1_mod,c(44,42,33))[[2]]
fmri_kimesurface(bigim1_mod,c(44,42,33))[[3]]
fmri_kimesurface(bigim1_mod,c(44,42,33))[[4]]

4 Interactive plotly Example

4.1 Plotly method: interactive way

The function fmri_image is used to create images for front view, side view, and top view of the fMRI image.

bigim1_mask<-bigim1_mod
for (i in 1:160) {
  bigim1_mask[,,,i]<-bigim1_mask[,,,i]*mask
}

Manually examples:

The plot without mask restriction.

fmri_image(bigim1_mod, option = "manually", voxel_location = c(40,40,30), time = 4)

The plot with mask restriction.

fmri_image(bigim1_mask, option = "manually", voxel_location = c(40,40,30), time = 4)

4.2 Forecasting with time series

The function fmri_ts_forecast will fit the ARIMA model for each location. This function is based on the TSplot_gen function from package TSplotly.

fmri_ts_forecast(smoothmod,c(41,44,33))

5 Motor area detection

5.1 fMRI data simulation

The function fmri_simulate_func is used to simulate the real-valued fMRI data with the specified dimension.

fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
                                   ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                                   dur = c(10, 10, 10, 10, 10, 10, 10, 10))

5.2 Stimulus detection

An integrated function fmri_stimulus_detect which can apply multiple methods is used to detect motor area.

“t-test” and “wilcoxon-test” can be applied to all real-valued fMRI data. “Hotelling’s T2”,“Wilk’s Lambda” and “gLRT” methods can be applied to all complex-valued fMRI data. “on_off_diff” and “HRF” methods can be applied to 4D real-valued fMRI data where the first method is calculating on-off difference period polar volume to get p-values and the second method is using hemodynamic response function and linear regression. The post-hoc stat-mapping filtering can be also applied to the calculated p-values if specified in the parameter or use the function fmri_post_hoc to do the process.

5.2.1 Examples

# statistical method HRF needs parameter ons and dur
pval1 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "HRF" , 
                             ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
                             dur = c(10, 10, 10, 10, 10, 10, 10, 10) )

# statistical method t-test for real-valued fMRI data
pval2 = fmri_stimulus_detect(fmridata= bigim1_mod, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))],
                             method = "t-test")

# statistical method Wilk's Lambda for complex-valued data
pval3 = fmri_stimulus_detect(fmridata = bigim1, mask = mask,
                             stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10)) ], 
                             method = "Wilks-Lambda" )

# do the fdr correction and the spatial clustering
# pval4 is the pval1 after the post-hoc processing
pval4 = fmri_post_hoc(pval1, fdr_corr = "fdr",
                                             spatial_cluster.thr = 0.05,
                                             spatial_cluster.size = 5, 
                                             show_comparison = FALSE)
summary(pval1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.8987  1.0000  0.8463  1.0000  1.0000
summary(pval2)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000096 0.8598920 1.0000000 0.8579245 1.0000000 1.0000000
summary(pval3)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000025 0.1877401 1.0000000 0.7476717 1.0000000 1.0000000
summary(pval4)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  1.0000  1.0000  0.9971  1.0000  1.0000

6 Motor area visualization

6.1 Visualization and comparison of p-value

6.1.1 3D visualization for p-value

pval1_3d = fmri_3dvisual(pval1, mask, p_threshold = 0.05, method="scale_p", multi_pranges=TRUE, title="HRF method")

pval1_3d$plot
pval4_3D = fmri_3dvisual(pval4, mask, p_threshold = 0.05, method="scale_p",
                         multi_pranges=TRUE, title="HRF method after post-hoc")

pval4_3D$plot

6.1.2 2D visualization for p-value

Generate the 2D plot of the p-value from sagittal, coronal and axial view.

Plot without hemodynamic contour.

for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=NULL, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}

Plot with hemodynamic contour.

for(axis in c("x", "y", "z")){
  axis_i = switch(axis, 
                  "x" = {35},
                  "y" = {30},
                  "z" = {22})
  print(fmri_2dvisual(pval1, list(axis, axis_i), 
                      hemody_data=hemody, mask=mask, 
                      p_threshold = 0.05, legend_show = TRUE, 
                      method = "scale_p",
                      color_pal = "YlOrRd", multi_pranges=TRUE))
}

6.2 Comparison of performance of different methods on the same fMRI data

6.2.1 3d p-value comparison

fmri_pval_comparison_3d is to visualize two p-value data to see how they perform differently on detecting stimulated parts by 3D plot. Here we compare the difference of stimulated parts of two different fMRI data with the same mask, since our original fMRI is too big to use here for different statistical tests.

fmri_pval_comparison_3d(list(pval1, pval2), mask, 
                                list(0.05, 0.05), list("scale_p", "scale_p"), 
                                multi_pranges=FALSE)

6.2.2 2d p-value comparison

fmri_pval_comparison_2d is to visualize whatever number of p values (generated by different statistical tests on the same fMRI data) to see their difference by 2D plot. For simplicity here we only compare two different 3D array p value data.

fmri_pval_comparison_2d(list(pval2, pval1), 
                               list('t-test', 'HRF'),
                               list(list(35, 33, 22), list(40, 26, 33)), 
                               hemody_data = NULL, 
                               mask = mask, p_threshold = 0.05, 
                               legend_show = FALSE, method = 'scale_p',
                               color_pal = "YlOrRd", multi_pranges=FALSE)

7 Three-phase ROI Analysis

7.1 Phase1: Detect Potential Activated ROI

phase1_pval = fmri_ROI_phase1(bigim1_mod, mask_label, mask_dict,
                              stimulus_idx = c(1:160)[rep(c(TRUE,FALSE), c(10,10))])$all_ROI$pval_t

7.2 Phase2: ROI-Based Tensor-on-Tensor Regression

phase2_pval = fmri_ROI_phase2(fmridata = bigim1_mod, label_mask = mask_label, label_dict = mask_dict, 
                              stimulus_idx = c(1, 21, 41, 61, 81, 101, 121, 141),
                              stimulus_dur = c(10, 10, 10, 10, 10, 10, 10, 10),
                              rrr_rank = 3, fmri.design_order = 2,
                              fmri.stimulus_TR = 3, method = "t_test")

7.3 Phase3: FDR Correction and Spatial Clustering

phase3_pval = fmri_post_hoc(phase2_pval , fdr_corr = "fdr",
                            spatial_cluster.thr = 0.05,
                            spatial_cluster.size = 5, 
                            show_comparison = FALSE)

7.4 3D visualization based on the activated areas by regions

fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, title = "phase1 p-values", rank = "value")
fmri_3dvisual_region(TCIU::phase1_pval, label_mask, label_index,
                     label_name, 5, title = "phase1 top five p-values", rank = "value")
# for 3D visualization, user needs to include empty region in the label
label_index = c(0, label_index)
label_name = c("empty", label_name)
fmri_3dvisual_region(TCIU::phase2_pval, label_mask, label_index,
                    label_name, title = "phase2 p-values")
fmri_3dvisual_region(list(TCIU::phase2_pval,TCIU::phase3_pval), label_mask, label_index,
                    label_name, title = "phase 2 & 3 p-values")
SOCR Resource Visitor number Web Analytics SOCR Email