Tutorial

1. Download and Install
 
Go to the download page and get the latest version of the PhLEM toolbox.  Extract the toolbox to the appropriate Matlab directory on your computer (e.g. C:\Program Files\MATLAB\R2007a) and add to your Matlab path.  

In the extracted PhLEM folder, there should be a directory called "demo".  In this directory you will find example datasets of ECG (.ecg), respiration (.resp), and pulse-oximetry/PPO (.puls) data taken off of a Siemens Trial 3T system. To work this tutorial, go into the demo directory and refer to these files.  

This program uses a lot of loader functions in the SPM5/SPM8 toolbox (with some attempt at backward compatibility to SPM2), so be sure to have the latest version of SPM installed and in your Matlab path.

2. Preprocessing Steps

See the README.txt file in the demo directory that comes included in the PhLEM.zip file.  This will walk you through some builds of simple noise models reported in the accompanying NeuroImage paper.  

Setting up the PhLEM Structure
You should be prompted to load each data file respectively.  This routine loads each dataset into Matlab, identifies the relative event peaks after smoothing/filtering, and transforms the inter-event interval into phase-time.  The loader functions are defined as strings in the PhLEM_defaults.m file.  To change which loaders you use, change the handle strings there.  Right now the loaders take a few minutes for datasets from typical experiments (e.g. 4-8 minute blocks), so be patient with the duration of this process.

Alternative Data Loader Functions
Note that if you do not have the physio log recording time-locked to your scan, I recommend using the Exvolt code developed by Geoff Aguirre and colleagues.  This does a beautiful job of time-locking data to your scan.  If you use this option, change the loader function handles in the PhLEM_defaults.m file to use the accompanying exvolt_data_loader.m function.

Filtering Options.
In terms of filtering, for ECG and PPG data, the toolbox uses a Butterworth bandpass filter to remove noise.  For respiration data, the toolbox uses a simple Gaussian kernel to low-pass filter.  Now, the default is to extract the high-frequency and low-frequency information from the PPG data separately.  Currently, this provides redundant phase information as the ECG and respiration signals respectively.  More on the relationship between signals will be discussed later. The default filtering parameters are all determined in the PhLEM_ana_opts.m file.  In some cases, you may require changing the band-pass filter window or guassian smoothing kernel depending on your particular datasets.  You can also point to your own parameters file (the demo includes a parameter file used in the published paper). .

The output is saved as a .mat file called 'test.mat'.  In most cases, this file would summarize all the phase information for the physiological signals.  To skip a data file, for example to ignore .ecg and .resp files to instead just use the PPO files, then just give only the file names of the data you are interested in  

>>PhLEM_setup('test.mat','ppgfile','FILENAME.PULS','ecgfile','FILENAME.ecg'); 

Calling the script this way will only process and output data for the pulse-oximetry and the ECG. As mentioned above, the output structure for the pulse-oximetry should have two fields.  One for the high-frequency PPO phase information and one for the low-frequency information.  The new version of PhLEM also includes an unfiltered PPG field that can be used to generate the down-sampled pulse-ox regressor model.

The best way to understand how this pre-processing works is to play around with the PhLEM_setup script with your own data and modify it to work properly for your datasets.  Make copies of it to setup your data as you wish or use the default version for quick, easy data analysis.

3. Generating Model Regressors

The routine that takes the high frequency phase information and puts it into usable form for the GLM is make_physio_regressors.m.  This function requires a string describing which signal to look at (e.g.- 'RESP','ECG','PPOHIGH', or 'PPOLOW') and the PhLEM data structure that is returned from the setup script.  It outputs a NxR array of values ("C")  where N = number of volumes per block and R = number of fourier expansions per data signal. In addition, the function outputs a cell array of regressor names ("names") for each column in C.

>>[C names] = make_physio_regressors('ECG',PhLEM);

If you are using the RETROICOR algorithm, then the C array is the average phase of each fourier series for the relevant data signal.  If you are using the 'rvhr' analysis option, this will generate a variation model.  If you are using 'raw' it will simply downsample your target signal.

The default TR for this function is 2 seconds.  To change this, provide the TR for your experiment as an optional input.  For example, if you had a 1 second TR, then call the function this way.

>>[C names] = make_physio_regressors('ECG',PhLEM,'TR',1);

The TR is always described in seconds but doesn't necessarily need to be in whole numbers (e.g. you can have a 1.5 sec TR).   See the README.txt file in the demo directory for more information.

4. Incorporating Phase Information Into Your GLM


The format of the output from make_physio_regressors is optimized for inserting into the SPM structure prior to running the modeling step in spm_spm.m.  To do this, simply add the C and names arrays for each block into the session field in the SPM structure. 

SPM.Sess(block_number).C.C = C;
SPM.Sess(block_number).C.names = names;

This will incorporate the physiological signals into your design matrix and treat them as covariates in the model.   The most optimal way to do this is to include this stage in any batch scripts you have running