Event-related responses in pupil size and fMRI measurements


It’s essential to understand the signals you’re looking at. The standard approaches to fMRI analysis assume the haemodynamic impulse-response function has a specific shape. Although slight changes in its timing and shape can be accounted for, it is nonetheless important to inspect the shapes of the event-related responses in one’s data. One way to correctly fit these event-related responses is to use linear regression to estimate the finite impulse-response caused by specific types of events.

Finite impulse response fitting has some advantages over and above event-related averages. If the impulse-response your data were filtered by is slow, i.e. a low-pass filter, it may be the case that responses to multiple events overlap. This causes problems for epoch-based event-related averaging, since some parts of the data will end up in the event-related responses of several different event types. Provided enough events of all types are present in the experiment (aided even more when inter-event periods are randomly drawn from an exponential distribution), deconvolution can take into account this overlap and correctly attribute the signal to specific types of events.

We have created a python class package, fir, to do this. This package makes it very easy to perform quite sophisticated analyses, distilling event-related responses from timeseries signals. See below for a brief example, for more elaborate examples, see this IPython notebook on GitHub

The package can be installed by issuing pip install fir --pre, where the --pre allows you to install software that’s still in development – as fir is.



from fir import FIRDeconvolution

# input_data is a numpy array sampled at signal_sample_frequency, for each of the separate signals (such as fMRI voxels). 
#     type: numpy array, (nr_signals x nr_samples)
# events are numpy arrays with event onsets in seconds. 
#     type: list of numpy arrays, (nr_event_types x nr_events_per_type)
# sample_frequency (Hz) is the frequency at which the original data are sampled (for example, 0.5 Hz for fMRI)
# deconvolution_frequency (Hz) is the frequency at which to perform the fit of the event-related signals
# the data will be resampled to this frequency before estimation.
# deconvolution_interval is a tuple, the interval in seconds over which to estimate the curves.

fd = FIRDeconvolution(
            signal = input_data, 								
            events = [events_1, events_2],						
            event_names = ['event_1', 'event_2'], 			
            sample_frequency = signal_sample_frequency,		
            deconvolution_frequency = deconv_sample_frequency,	
            deconvolution_interval = (-5,15)		
            )

# we then tell it to create its design matrix
fd.create_design_matrix()

# perform the actual regression, in this case with the standard numpy.linalg backend
fd.regress(method = 'lstsq')

# and partition the resulting betas according to the different event types
fd.betas_for_events()

fd.calculate_rsq()


And this code might return the type of curves shown in the above figure. We encourage you to give this package a try when you have fMRI data from voxels, ROIs, or pupil size timeseries data.