HRfunc

Getting Started



The HRfunc tool can be used to estimate hemodynamic response functions and neural activity in functional near infrared spectroscopy (fNIRS). Estimate neural activity in each subject and a group level hemodynamic response function (HRF) through Python and MNE following these steps. Check the Q & A or reach out if you have any questions.

  1. Prepare raw fNIRS scans and events
  2. (Optional) Estimate subject-level HRF estimates for each subject
  3. (Optional) Generate a subject pool HRF distribution
  4. Estimate channel-wise neural activity for each subject


Video guide of how to use HRfunc to estimate HRFs and neural activity...



Step 1 - To first estimate subject level HRFs, load in your fNIRS data through MNE. Initialize a hrfunc.montage() for estimating neural activity, alongside HRFs if you have an event-based experiment design. 

import mne
import hrfunc as hrf

# Load in NIRX scans
raw_scan = mne.io.read_raw_nirx("path/to/your/nirx/scan/")

# Load a hrfunc montage
montage = hrf.montage()


Step 2 (Optional) - To then estimate your own HRFs, if you have an event timeseries tied to your scans, load in your events (a list of 0's and 1's indicating when events occured in your fNIRS scan) however you store them. In this example we stored our events in a .txt file. Using these events alongside their relative fNIRS scans, iterate through each of your raw scans passing them into the montage.estimate_hrf() function. The HRfunc tool will automatically preprocess them for deconvolution, which is distinctly different than standard hemoglobin preprocessing. The lmbda (lambda) variable can be tuned to increase and decrease noise suppresion within the signal.

NOTE: If you did preprocess your data for deconvolution already, set the parameter preprocessed = True and HRfunc will skip preprocessing before estimating HRFs.


# Load in events and format as a list of 0's and 1's representing an event sequence in NIRX sampling space
with open("path/to/events.txt", "r") as file: # Or however you save/load events
    events = [int(line.split('/n')[0]) for line in file.readlines()]    

# Iterate through each of your MNE NIRX scans and estimate an HRF for it
for scan in raw_scans:
    montage.estimate_hrf(scan, events, duration = 30.0, lmbda = 1.0)


Step 3 (Optional) - Once each scan has been passed into your montage to have HRFs estimated, you can now average across subject estimates to generate a subject-pool wide HRF estimate for each channel. After generating a subject-pool wide HRF distribution for each NIRX channel, an HRF estimate will be attached to each of your hrf montage channel nodes for use in estimating neural activity. Save the object as a json object to use with other study analysis, share with collaborators, or share with the wider neuroimaging community. 

# Generate subject pool wide distribution
montage.generate_distribution() 
montage.save('study_HRFs.json') 

# Load previously estimated HRF montage back up for step 4 when ready
montage = hrf.load_montage('study_HRFs.json')


Step 4 - Now using our estimated HRFs we can estimate neural activity! To accomplish this you must iterate through each of your scans and pass them into your montage.estimate_activity(raw). This will replace channel hemoglobin data in place using using MNEs raw.apply_function() to each channel using you're estimated channel HRF.

NOTE: You do not need to estimated HRF's to estimate neural activity. Without using estimated HRF's, a canonical double-gamma HRF is used for deconvolving each channel. Accuracy of HRFs may vary depending on how well the HRF matches the true HRF for your subject pool. 

# Estimate neural activity in place and save the nirx object
for scan in raw_scans:
    scan = montage.estimate_activity(scan)
    scan.save('sub-123_deconvolved.fif')