Physiological Artifact Removal Tool (PART)
Physiological cycles such as breathing and heart rate can add noise to 4-dimensional MRI data. For example, voxels in a fMRI acquistion may become brighter or darker depending on the cardiac phase. The statistics used to analyze fMRI data is based on the ratio of predicted signal versus noise, so if we can reduce the noise in our data, we can effectively increase our ability to detect brain activation. At first glance, one would think that you could simply apply a ‘no’tch filter to eliminate the frequency of the noise we are interested in. For example, the fMRI signal has a cycle of ~0.04 Hz (6 sec to peak, 12 sec to baseline, 18 sec maximum undershoot and 24 sec back to baseline) while the heart rate is about 1 Hz and the respiration cycle is about 0.3 Hz. Indeed, for block fMRI designs, I recommend applying low pass temporal filtering to your data (this option is on the ‘pre-stats’ tab of FSL’s FEAT). However, there are two reasons why this solution is less than ideal. First, with event related fMRI designs, some of the signal may be at higher frequencies. Second, we have the problem of aliasing – if a frequency is faster than half our sampling rate, it can appear as a slower frequency (the Nyquist theorem). Aliasing is the reason that wheels can appear to move backwards when we watch TV: the TVs 60Hz sampling rate is too slow for the wheel’s angular rotation. For example, consider if we sample a signal at 1.5 times it’s period. In the graph below the black line shows the actual signal fluctuation, while the red dots show our samples. Observing our samples, we would infer that the signal was half it’s actual frequency (dotted line):
Note that with fMRI we typically observe a voxel of the brain around every 2-3 seconds, slower than the critical Nyquist frequency for cardiac and respiration. If we sample every 2 seconds, we can not observe frequencies with a period shorter than 4 seconds, while the heart rate is approximately once every second. In short, physiological artifacts can add noise that reduces our ability to detect brain activation. This problem is especially acute for resting state fMRI, where aliased physiological signals could mimic real activation (while with standard fMRI we are phase locking the signal to external events, so physiological signals may reduce our ability to detect signals but should not cause false alarms). There are a couple of solutions that can help reduce the effect of physiological noise:
- Measure physiological cycles during the MRI scan, and then retrospectively remove changes in image intensity that are phase-locked to the physiological cycle. Note that the physiological variables can either be covaried out during analysis (Corfield et al. 1999, J Appl Physiol. 86: 1468-1477) or by preprocessing the data. Examples of this latter method have been described by Glover et al. 2000 Magn Reson Med. 44:162-167 and Deckers 2006 Neuroimage. 233: 1072-1081. Perhaps the most popular algorithm is Glover’s RETROICOR, which is inlcuded as a plugin for AFNI.
- Acquire volumes faster than the nyquist frequency, in other words have a TR <0.5sec. For example Anand et al.(2005) Neuropsychopharmacology 30, 1334-1344 use a TR of just 0.4 sec, collecting just four slices. This approach does include drawbacks – the short TR means there will be very little signal to measure (as there is little T1 recovery), furthermore, the small number of slices acquired necessarily mean you will only see a small part of the brain (and you can only observe effects in parts of the brain that you measure), furthermore, the small number of slices means you can not normalize your images to conduct conventional group analyses. Therefore, while this is an excellent way to model and eliminate physiological noise, it requires careful analysis – the paper by Anand et al. demonstrates how careful region of interest analysis can effectively use this technique.
- The whole brain’s fMRI signal fluctuates with physiological (respiratory) cycle. Therefore, one approach is to model this effect as a regressor in your analysis (Birn et al. 2006 Neuroimage. 31:1536-1548; Vincent et al. 2006 J Neurophysiol. 96:3517-3531). This might be a useful technique for resting state fMRI. However, for event-related fMRI or situations where large portions of the brain cycle together, I wonder if this approach might have some of the drawbacks as global scaling.
- Acquire data with the acquisition gated to the cardiac cycle. Unfortunately, the TR (time between volumes) would vary between volumes, so the signal (T1-effects) would also change. You would be replacing one type of noise with another. In theory, you could model the known modulation in TR. I do not know of any groups that have used this approach.
- While conventional fMRI scans take around 2-4 seconds per volume, the amount of time to collect each slice is around 50-100ms. In theory, it should be possible to use the 3D information to both infer and remove physiological artifacts. This would be a truly elegant solution, as it would not even require physiological monitoring. While this option has been discussed in theory, I do not know of any implementations.
Introspectively, one might think that heart rate has a larger influence on fMRI signal than respiration (as fMRI measures changes in blood flow). However, the reverse is true: heart rate only has localized effects due to the blood vessel’s pulsatility (Dagli et al 1999 Neuroimage 9: 407-415), while respiration tends to correlate with head motion that influence the magnet’s homogeneity typically resulting in global effects.
I have adapted a simple technique described by Deckers (see above) to reduce physiological noise. This alogorithm divides the period between successive pulses/breaths into a number of bins. The observed data is then classified according to which bin it is closest to. The mean signal for each bin is computed, and this observed fluctuation is then removed from the data. Following Deckers, my software uses 40 bins for cardiac data and 20 bins for respiration bins. Deckers. In the algorithm specificied by Deckers, the bins span the mean physiological period plus two standard deviations (in other words, if the mean time between cardiac pulses is 1.02sec +/-0.1sec, then the bins span 1.22sec). In my algorithm I use the Median plus 1.5 inter-quartile ranges. I did this as participant movements in the scanner can lead to sensors failing to detect a few heart beats – and the mean and standard deviation are influenced by these missing samples. Furthermore, my software will interpolate up to two missing heart beats.
To use my software you need to record the pulse and/or the respiration. This is easy to do with modern Siemens scanners, which can generate recording files named *.puls and *.resp (Philips scanners can collect .log files). With other scanners, you can use my software if you generate FSL-format 3 column text files where the first item in each column reports the physiological trigger time (e.g. 0.12 1 1 1.21 1 1… if the first heart beat occurred 0.12 seconds after MRI acquisition and the second pulse occurred 1.21 seconds into the recording). Below is a sample of my own breathing (blue line) and heart rate (red line), with the triggers shown in magenta (heart) and cyan (breating).
Here is what you will need to run my software, you can get a sample dataset (23mb) by clicking here (this also includes the optional physiological data viewer software shown above, which requires the Windows operating system):
- Your physiological data (either Siemens .puls/.resp, Philips .log or FSL 3-column text files).
- Your MRI data (a 4D NIfTI image).
- My software currently assumes ascending acquistion with continuous sampling (e.g. do not use this with sparse fMRI).
- Your data should be motion corrected prior to physiological correction (e.g. use FSL’s mcflirt).
- A copy of NPM.
Simply Launch NPM and select Utilities/PhysiologicalCorrection.. The choose your physiological and MRI data. The processing requires about a minute, with a new image generated with the prefix ‘i’, so if you processed the image c:\temp\rchris.nii.gz you would get an image named c:\temp\irchris.nii.gz. The software generates information that tells you how many samples were in each bin. The image below shows a fourier transform for a cortical boxel before (black) and after (red) correction – note that this voxel shows little change in signal.
A sample output will look like this:
Physiological Artifact Removal Tool started = 10:02:53 AM
Assuming continuous fMRI ascending acquisition with TR = 2.1000sec
Correction for C:\cygwin\home\mscae\20061220_140508\physiolog_WED_DEC_20_14_39_34.resp
Time per vol (TR) [sec] 2.1000
Triggers n/First…Last [vol] 229.00/0.05..346.72
Q1/Median/Q2 [sec] 2.92/3.16/3.44
Bin n/Range [sec] 20/-2.36…2.36
voxels without variance (outside brain) %: 76.36
voxels with variance which were corrected %: 98.77
Bin 1 -1.01 0
Bin 2 -0.89 2
Bin 3 -0.77 1
Bin 4 -0.65 10
Bin 5 -0.53 19
Bin 6 -0.41 25
Bin 7 -0.30 32
Bin 8 -0.18 14
Bin 9 -0.06 33
Bin 10 0.06 31
Bin 11 0.18 25
Bin 12 0.30 24
Bin 13 0.41 30
Bin 14 0.53 29
Bin 15 0.65 24
Bin 16 0.77 22
Bin 17 0.89 17
Bin 18 1.01 7
Bin 19 1.12 0
Bin 20 1.24 0
Input = C:\cygwin\home\mscae\20061220_140508\rachris.nii.gz
Output = C:\cygwin\home\mscae\20061220_140508\irachris.hdr
Physiological Artifact Removal Tool finished = 10:03:16 AM
Ensure that your fMRI protocol has been adapted to save physiological data. On the Siemens Trio TIM system, use the wireless finger clip to record pulse data – the laser light should illuminate the proximal portion of the finger nail (not the finger pad, as claimed in the manual), while you can record respiration with the respiration cushion (use the velcro chest strap to hold the cushion snug, but loose enough that the cushion expands and contracts with each breath). You can set up the Siemens console to show the pulse and respiration simultaneously (assuming you are using System B13 or later), with the resulting log files saved in c:\medcom\log (you will need to activate advanced user mode to view, copy and delete these files).