1 Preface: Project outline

1.1 fMRI time-series analysis and Primary Motor Area (PMA) detection

Functional magnetic resonance imaging (fMRI) is one of the most widely used tools to study the neural underpinnings of human cognition. One of the most important application for fMRI time-series data is to detect the Primary Motor Area (PMA) of the brain. My goal in this report is to develop several simple and general methods to detect PMA based on a given brain scan experimental design (like the following finger tapping experiment in this report) and the corresponding fMRI data. A typical fMRI data is a 4 dimension time-series data consists of 3 spatial dimensions (X, Y, and Z axes) representing the brain area. Each X, Y, and Z coordinate represents a cuboid element (a brain area) of the fMRI scan result, we call each of them as a “voxel”. The last dimenson is a temporal dimension representing the time when the MRI machine conducts a brain scan on every voxel, which makes each voxel contain a time-series data which mainly captures the changes in blood flow and oxygenation associated with neural activity (the hemodynamic response), and on the differing magnetic properties of oxygenated and deoxygenated blood of that certain voxel. In particular, we call the signal captured by the time-series data in each voxel as the blood oxygen level-dependent (BOLD) signal.

For any certain given biain scan experimental design, we would ask our subjects to conduct that certain experiment with a certain time, we call all MRI scans conducted in those periods as “ON epochs”. By contrast, there are also certain time that our subjects are asked to rest without any certain movement. Those periods would act as the control group with the “ON epochs”, which we call as “OFF epochs”. Moreover, for a given experiment, we would expect that only a few certain brain areas would be activated, especially for a given experiment related to the movement of a certain human body, we would expect activation of one or several certain areas in the Primary Motor Area (PMA) of the brain.

Primary Motor Area (PMA) of the brain

Primary Motor Area (PMA) of the brain

Accurately detecting the Primary Motor Area of a certain experiment would be helpful for neural science researchers to research on the functionality of a human brain and neural activities. In this project, I try to achieve this goal by applying time series models on the time series data of each voxels. A straightforward intuitive idea is that for voxels in PMA, the BOLD signals for ON epochs and OFF epochs would be different. Thus would show some particular or seasonal pattern. (BOLD signal is high for ON epochs and low for OFF epochs). And voxels that are not activated wouldn’t show a clear pattern (may be white noise signal). And I will try to develop a detecting mechanism based on this idea.

1.2 Data description

In this report, the data I use comes from a sensitive source. It is stored in the file “MOD1.RData”. It contains the fMRI data which comes from the following finger tapping experiment:

Finger-Tapping Experiment

Finger-Tapping Experiment

My data is a 4D fMRI data with \(64*64*40\) voxel locations (i.e. Brain scan area: X axis: 64, Y axis: 64, Z axis: 40). Each voxel location has 160 scans arranged by time. There are 3 seconds pause between each scan. Thus each experiment would last for 8 minutes. Just like what is illustrated in the graph above, during those 160 scans there are alternant ON and OFF epochs, each ON and OFF epochs last for 30 seconds, thus give us 8 ON and OFF alternant periods for the whole experiment.

All analysis in this report is based on this fMRI data. But this report aims to introduce a simple and general method that could be applied in any similiar MRI scan experiment. More illustrations will be shown in the following report.

2 Data visualization

The following graph shows the time series data of two voxels. First voxel is located at X=44, Y=40, Z=32. And the other one is located at X=10, Y=10, and Z=10.

To visualize more detailed information of each voxel, a function “plotfMRI_TS” is developed. By key in the spatial and temporal information of a certain voxel, researchers would see the time series data and the voxel location on each axis. Moreover, two smoothing trends are also added in the time series plot, they are both achieved by loess smoothing. One with span=0.5, intends to show the trend of this time series. The other with span=0.05, intends to show the noise.

Below is the visualization result of voxel X=41, Y=40, and Z=32. In the upper left time series plot, blue line is the raw time-series data, red line is for the noise signal, and the green line is for the trend signal. The upper right plot is the top view of the brain (Z axis view), the left bottom plot is the side view (X axis view) , and the right bottom is the front view (Y axis view).

This particular voxel location is chosen to show the time series signal of a typical activated location a brain in this finger tapping experiment. This will serve as a typical example in the following report.

By contrast, the following graph comes from voxel location X=10, Y=10, and Z=10. It is chosen to show the typical signal of an unactivated area in this MRI experiment. In fact, based on the visualization result, you can easily see that this is a voxel area outside the brain. Still, there are some BOLD signal detected in those areas. It is most likely caused by the noise when conducing a MRI scan. When comparing this data with the data inside the brain area, we can easily see that the intensity of the BOLD signal outside the brain area is much smaller than signal inside the brain area. This could serve as one method to rule out unactivated voxel locations outside the brain area. However, in this report, I will make the judgement based on the trend of each voxel time series, as this method is more general and could be applied to voxels inside the brain area.

This particular voxel location, X=10, Y=10, and Z=10 will also serve as a typical example in the following report.

3 Time series analysis

3.1 domain analysis

As our finger tapping experiment is based on 8 periods of ON and OFF epochs, each period contains 20 brain scans. Our first aim is to detect that whether a certain voxel would present a clear periodical behavior.

## [1] "Dominant period for voxel X = 41 Y = 40 Z = 33 is 22.8571428571429"

Voxel X=41, Y=40, and Z=33 shows a clear seasonal trend. Based on the smoothed periodogram, the dominant period given is around 23, which is close to the expected period (20) in this experiment.

## [1] "Dominant period for voxel X = 10 Y = 10 Z = 10 is 4.57142857142857"

By contrast, the periodogram given by the unactivated voxel shows that there is not a very clear seasonal trend for unactivated voxels. And the detected dominant period doesn’t make much sense.

Based on the results above, we could apply the periodogram graph as a intuitive way to detect PMA.

3.2 Smoothing: filter trend, noise, and cycle

By loess smoothing, we could decompose the trend, noise, and cycle signal for each voxel.

Based on the decomposition of voxel X=41, Y=40, and Z=32, we can see that we might detect a clear linear trend of the time-series data of each voxel. But it is hard to detect any seasonal trend based on the ON and OFF epochs by this decomposition. Especially, the noise signal doesn’t really make sence as it is too close to the original time-series data. This might be caused by the reason that we only have 160 time points in each voxel, and it is hard to detect local noise based on such a small sample size. Thus the cycle signla also doesn’t make much sense here.

As for the unactivated voxel, the voxel and cycle signal also don’t make much sense. And here we could find that there is not a very clear linear trend in this voxel. This is reasonable as we might expect that there is only white noise signal in unactivated voxel. However, we shouldn’t be too suprised if we detect a linear trend for those voxels.

Based on the analysis above, for fMRI data with small sample size in each voxel, it is hard to detect the cycle signal based on the loess smoothing decomposition method. But we could still detect the linear trend based on this method. This also suggests that to truely if a voxel is an unactivated voxel, we might wish to detrend first and then judge that if the rest signal is stationary or even just pure white noise.

3.3 Trend analysis: linear, quardatic, and cubic trend

The following part is to analysis that whether we might think about introducing higher order trends (quadratic, cubic trend) inside our time series model.

The ARMA model is selected based on the AIC table. And all three trends information is added to the ARMA model as external regressors.

MA0 MA1 MA2 MA3 MA4
AR0 1558.825 1558.527 1555.790 1556.679 1555.139
AR1 1557.712 1557.907 1553.335 1554.910 1557.131
AR2 1555.757 1552.795 1554.795 1556.698 1558.674
AR3 1556.319 1554.795 1556.792 1535.388 1558.445
AR4 1555.539 1556.515 1558.515 1558.225 1560.134 
## 
## Call:
## arima(x = bigim1_mod[41, 40, 32, ], order = c(2, 0, 1), xreg = X_trend)
## 
## Coefficients:
##           ar1     ar2     ma1  intercept  linear   square  cubic
##       -0.6057  0.2610  0.7694  1275.6156  1.4492  -0.0169  1e-04
## s.e.   0.1467  0.0847  0.1428    12.4508  0.6795   0.0100  1e-04
## 
## sigma^2 estimated as 867.6:  log likelihood = -768.4,  aic = 1552.79
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE     MAPE      MASE
## Training set 0.1232036 29.45543 23.47029 -0.04110395 1.794745 0.7470589
##                     ACF1
## Training set -0.00493342

Based on the model, we see that only the linear trend is significant.

Based on the graphs given above, we see that this ARMA model with three trends could caputure the main characteristics of the activated voxel. The residuals left seem follow the normality assumption and doesn’t show a clear lag. However, there still seems to have a dominant period similiar to the original data left in residuals, and the quadratic trend and cubic trend are not significant. This suggests that we might wish to change another common trend to represent the activated voxels.

MA0 MA1 MA2 MA3 MA4
AR0 1423.481 1425.452 1425.748 1426.735 1428.301
AR1 1425.458 1418.305 1419.546 1428.676 1426.979
AR2 1426.004 1419.568 1420.149 1421.237 1419.994
AR3 1426.160 1421.513 1422.036 1413.975 1420.476
AR4 1427.817 1423.194 1423.073 1421.523 1414.943 
## 
## Call:
## arima(x = bigim1_mod[10, 10, 10, ], order = c(0, 0, 0), xreg = X_trend)
## 
## Coefficients:
##       intercept   linear  square  cubic
##         47.3073  -0.1030  -2e-04  0e+00
## s.e.     6.4917   0.3481   5e-03  1e-04
## 
## sigma^2 estimated as 402:  log likelihood = -706.74,  aic = 1423.48
## 
## Training set error measures:
##                        ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 2.359712e-12 20.04923 16.5211 -63.66982 88.64721 0.7188974
##                     ACF1
## Training set -0.01217592

For the unactivated voxel, all three trends are not significant, and we are suggested to apply an ARMA(0,0) model, which suggests that the original time series data might just be a white noise signal.

4 Time series models with BOLD response signal

4.1 Simulate BOLD response signal & BOLD signal regression model.

As we expect that for voxels inside the Primary Motor Area would show a different seasonal trend for ON and OFF epochs, we would expect that there is a common trend that is shared by all activated voxels that could represent the typical BOLD signal of activated PMA voxels. Based on package “fmri” and article “A General Statistical Analysis for fMRI Data”, we find out that R could easily generated this certain BOLD response signal based on functions “fmri.stimulus” and “fmri.design”. One more advantage we have here is that these two functions are also general methods that could apply to similiar experiments. We only need to know the particual experimental design to stimulate the expected BOLD response of activated voxels for each experiment.

The following graph shows the expected BOLD response of this finger tapping experiment, along with the linear signal. Here we keep the linear trend as analysis above shows that we might face a significant liear trend in each voxel. These two trends would act as the external regressors for the time series models

If we only combine those two signals above and fit a regression model with the time-series data, we would get the following result. We see that it could fit the overall trend well but it might not able to capture other time series information from the original time series data. Thus we might conside ARMA model with BOLD and linear trends.

5 General methods for detecting Primary Motor Area

Combining all methods we applied above, we could apply the following methods to detect whether a voxel is inside the Primary Motor Area or not.

The first method is a preliminary screening method aiming to detect voxels that couldn’t be inside the PMA.

  1. Plot the smoothed periodogram to detect whether the dominant period is close to the alternant ON and OFF epochs period.
  2. If each voxel would have sufficent time series samples, we could think about conducting trend, noise, and cycle decomposition to detect potential cycle trend.
  3. Fit an ARMA model with only linear trend by AIC table. An evidence for a voxel is not inside the PMA is that we get an ARMA(0,0) model with or without a significant linear trend.
  4. Fit a linear model with linear and BOLD response. An evidence for a voxel is not inside the PMA is that the BOLD signal coefficient is not significant in this linear model.

A function “fMRI_TS_overview” is developed to conduct all procedures above. And if we would go to the next procedure to detect a certain voxel in more detail, an AIC table would be given based on the ARMA model with both linear and BOLD response. And we could select the best ARMA model that would be applied in the next procedure.

For voxels that might be inside the PMA, a function “fMRI_TS_analysis” is developed to conduct more detections.

  1. The best ARMA model with both linear and BOLD signals are fitted. The dominant factor of detecting the PMA is to see whether the coefficient of the BOLD signal inside the ARMA model is significant. As if we get a significant BOLD signal, it represents that the time series data in this voxel could fit perfectly in the expected BOLD signal response, thus represent a high possibility of being an activated voxel.
  2. Diagnosis of residuals are conducted. If the residuals are not close to white noise, we might need to dig deeper into that voxel to determine whether it is inside the Primary Motor Area.

5.1 Examples and demos for the general Primary Motor Area detecting procedures.

5.1.1 Example 1: voxel X=32, Y=32, Z=20

This voxel is chose as it is the center of the fMRI data. This voxel must be inside the brain but according to the knowledge about Primary Motor Area, it is not one of the part of the PMA.

## [1] "Dominant period for voxel X = 32 Y = 32 Z = 20 is 53.3333333333333"

The dominant period doesn’t match the expected dominant trend.

MA0 MA1 MA2 MA3 MA4
AR0 1554.957 1555.841 1554.738 1554.493 1556.074
AR1 1555.547 1556.536 1554.256 1556.212 1557.228
AR2 1555.272 1554.727 1555.169 1556.314 1559.016
AR3 1555.083 1556.293 1556.173 1557.681 1559.939
AR4 1555.518 1556.946 1558.900 1560.672 1562.148

AIC table given by ARMA model with only linear trend suggests that after detrending, the fMRI signal is highly likely to be white noise.

## 
## Call:
## lm(formula = Y ~ ., data = voxel1_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -64.831 -22.289  -2.446  21.641  78.667 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1129.022      4.912 229.855  < 2e-16 ***
## BOLD_signal      8.434     12.000   0.703    0.483    
## linear_signal  -30.346      4.263  -7.118 3.68e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.86 on 157 degrees of freedom
## Multiple R-squared:  0.2494, Adjusted R-squared:  0.2398 
## F-statistic: 26.08 on 2 and 157 DF,  p-value: 1.666e-10

The BOLD signal is not significant in the linear model.

No interesting information we could get from the decomposition plot.

Based on all information we gathered above, we could conclude that voxel X=32, Y=32, and Z=20 is not inside the PMA.

5.1.2 Example 2: voxel X=41, Y=24, Z=32

This voxel is in the different Y axis position with the activated voxel X=41, Y=40, Z=32. It is in the back of the brain and in theory it is a typical unactivated area.

## [1] "Dominant period for voxel X = 41 Y = 24 Z = 32 is 3.26530612244898"

The dominant period doesn’t match the expected dominant period.

MA0 MA1 MA2 MA3 MA4
AR0 1504.934 1504.725 1506.672 1507.506 1509.481
AR1 1504.880 1506.705 1506.755 1508.284 1507.741
AR2 1506.435 1508.000 1510.106 1508.627 1503.615
AR3 1507.176 1509.173 1504.806 1504.168 1503.578
AR4 1509.170 1507.339 1509.041 1505.799 1507.629

AIC table given by ARMA model with only linear trend suggests that after detrending, the fMRI signal is highly likely to be white noise or could be fitted by ARMA(0,1) model. Actually, there is not a very strong evidence that this is not a white noise signal after detrending.

## 
## Call:
## lm(formula = Y ~ ., data = voxel1_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -58.859 -21.453   2.991  17.782  74.909 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1586.299      4.204 377.368  < 2e-16 ***
## BOLD_signal     -5.668     10.270  -0.552    0.582    
## linear_signal  -29.730      3.648  -8.149 1.09e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.41 on 157 degrees of freedom
## Multiple R-squared:  0.2972, Adjusted R-squared:  0.2883 
## F-statistic:  33.2 on 2 and 157 DF,  p-value: 9.439e-13

The BOLD signal is not significant in the linear model.

No interesting information we could get from the decomposition plot.

Based on all information we gathered above, we could conclude that voxel X=41, Y=24, and Z=32 is not inside the PMA.

5.1.3 Example 3: voxel X=42, Y=41, Z=33

This voxel is very close to the activaed voxel. And it is very likely to be a voxel inside the PMA as PMA should be a few concentrated area.

## [1] "Dominant period for voxel X = 42 Y = 41 Z = 33 is 20"

The dominant period matches the expected dominant period and the ACF plot shows a trend pattern.

MA0 MA1 MA2 MA3 MA4
AR0 1564.037 1561.601 1561.123 1554.182 1554.771
AR1 1560.357 1557.559 1557.398 1555.085 1555.664
AR2 1558.528 1558.348 1556.825 1556.477 1557.629
AR3 1555.289 1557.078 1554.858 1554.408 1560.469
AR4 1556.916 1556.201 1547.672 1549.740 1558.260

AIC table suggests an ARMA(0,4) model.

## 
## Call:
## lm(formula = Y ~ ., data = voxel1_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -88.951 -15.472  -0.669  19.967  67.789 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1682.656      4.617 364.477  < 2e-16 ***
## BOLD_signal     63.497     11.279   5.630 8.12e-08 ***
## linear_signal    6.284      4.007   1.568    0.119    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29 on 157 degrees of freedom
## Multiple R-squared:  0.1739, Adjusted R-squared:  0.1634 
## F-statistic: 16.53 on 2 and 157 DF,  p-value: 3.066e-07

Clear signifance on the BOLD signal coefficient.

All information we gathered above suggest that this is a possible voxel inside the Primary Motor Area. Thus more detection should be conducted.

MA0 MA1 MA2 MA3 MA4
AR0 1536.617 1538.284 1538.233 1533.819 1535.819
AR1 1538.212 1537.660 1537.304 1535.819 1537.722
AR2 1538.132 1538.147 1531.779 1537.354 1538.169
AR3 1534.845 1536.460 1537.942 1539.207 1541.195
AR4 1536.184 1537.954 1539.038 1541.116 1533.488

Using this AIC table, we get that the best ARMA model with both BOLD and linear trend is ARMA(2,2)

## 
## Call:
## arima(x = data_fmri, order = c(p, d, q), xreg = X_fmri)
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  intercept  BOLD_signal
##       1.3614  -0.9414  -1.3347  1.0000  1682.4905      59.0296
## s.e.  0.0392   0.0369   0.0300  0.0319     4.9634      10.9860
##       linear_signal
##              6.3915
## s.e.         4.3081
## 
## sigma^2 estimated as 740.9:  log likelihood = -757.89,  aic = 1531.78
## 
## Training set error measures:
##               ME RMSE MAE MPE MAPE
## Training set NaN  NaN NaN NaN  NaN

After applying the best ARMA model, coefficient of the BOLD signal (based on z-test) is significant. This is a good evidence to prove that this is a voxel inside the PMA.

Diagnosis on residuals show that the residuals are highly likely to be white noise. Thus prove that the ARMA model we developed above are reasonable.

Thus we could conclude that the voxel X=42, Y=41, and Z=33 is inside the PMA.

5.2 A bold approach to apply a simple detecting method to all voxels

As we could see that the significance of the coefficience of BOLD signal inside the ARMA model with both linear and BOLD signal is a leading feature to determine whether a voxel is inside the PMA or not. Thus we could conduct a simple z-test on all voxels to calculate the p-value of the coefficient of BOLD signal for all voxels. This is not a very appropriate approach. There are two major flaws for this method. Firstly is that we might detect too many activated voxels. To solve this, we might apply FDR correction method or Benjamini-Hochberg Procedure to control detected voxels. Secondly, the detected voxels are likely to be scattered and not concentrated on several major area. To solve this issue, we might think about applying spatial clustering to concentrate the PMA we detected. However, this is not the only issue we face. As we didn’t take the correlations between voxels into consideration. Some articles have addressed the inaccuracy for detecting the Primary Motor Area only based on single voxel. And we might be able to get a better method by some bayesian analysis approaches. But this is already beyond the discussion of this report. Nonetheless, this bold approach would serve as a good tool to have a overall visualization on the detected primary motor area.

The concentrated red areas are the detected Primary Motor Area. But clearly, more visualizations could be done based on the detected p-value. We would be able to show a very nice 3D visualization of the brain primary motor area after conducing FDR correction and Spatial Clustering.

6 Reference

  1. Article: Statistical Analysis of fMRI Time-Series: A Critical Review of the GLM Approach
  2. Article: A General Statistical Analysis for fMRI Data
  3. Article: Analysis of fMRI Timeseries: Linear Time-Invariant Models, Event-related fMRI and Optimal Experimental Design
  4. Article: Bayesian Models for fMRI Data Analysis
  5. Book: Spectral Analysis of fMRI Signal and Noise
  6. Slides: The general linear model for fMRI
  7. Previous project:Monthly Traffic Fatalities in Ontario (1960-1974)
  8. Class notes and previous exams
  9. Massive assistance from my friends Zhe Yin and Jinwen Cao
  10. Suggestions and guidance from my professor Ivo Dinov.