Cerebrovascular Reactivity Mapping with phys2cvr

Cerebrovascular Reactivity Mapping with `phys2cvr`#

Cerebrovascular Reactivity (CVR) Mapping can greatly benefit from multi-echo (ME) data acquisition and independent component analysis (ICA) denoising, especially when it is based on breath-holds, cued deep breaths, or resting-state.

While ME-ICA will greatly reduce motion-related noise (see Moia et al. [2021]), because CVR is an increase in blood flow triggered by vessels dilation, there is a component of signal that is proton-density related. For this reason, the default decision tree in tedana may categorise certain components that should (or at least could) be kept as noise, so it is strongly suggested to perform a manual revision of the automatic categorisation and/or adapt the decision tree to accept a better suited set of components.

On top of that, in the context of breath-holds and cued deep breaths, attention should be made in the denoising stage, as an aggressive approach will remove all sources of signal that inform CVR and the associated haemodynamic lag. Instead, a conservative denoising approach, where noise IC timeseries are orthogonalised to the (temporally shifted) \(P_{ET}CO_2\) and the non-noise IC timeseries is to be preferred (Moia et al. [2021]).

Assuming that you are following the other tutorials in this notebook, once you obtained your optimally combined signal and your ICA decomposition that has been manually inspected (and corrected if needed), you will need to make sure you have:

a whole brain mask in functional space that will confine phys2cvr operations within it - the desc-adaptiveGoodSignal_mask.nii.gz output of tedana can work fine, but if you have a different mask that is also ok;
if you want, a ROI mask in functional space which average signal will be used as reference to align the regressor (e.g. the \(CO_2\) trace) and which average haemodynamic lag will be set as 0-lag - grey matter or cerebellum are commonly used as reference ROI;
if you want to run ICA denoising, the timeseries of rejected and accepted components from tedana in two different files; the easiest way to obtain these are using the mixing matrix from tedana and the manual classification from rica with the help of pandas:

import pandas as pd

# Assuming RICA was used for manual classification, read the downloaded file (adjusting path as necessary)
man_class = pd.read_csv('../manual_classification.tsv', sep='\t', header=0)

# Read the ICA mixing matrix  (adjusting path as necessary)
ica_mix = pd.read_csv('../desc-ICA_mixing.tsv', sep='\t', header=0)

# Extract the list of rejected and accepted components
rej_comp =  man_class[man_class['classification'] == 'rejected']['Component'].tolist()
acc_comp =  man_class[man_class['classification'] == 'accepted']['Component'].tolist()

# Extract rejected vs accepted timeseries and save them
rej_ts = ica_mix[rej_comp]
acc_ts = ica_mix[acc_comp]

acc_ts.to_csv('../accecpted_ic_timeseries.csv', index=False, header=False)
rej_ts.to_csv('../rejected_ic_timeseries.csv', index=False, header=False)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[1], line 4
import pandas as pd

# Assuming RICA was used for manual classification, read the downloaded file (adjusting path as necessary)
----> 4 man_class = pd.read_csv('../manual_classification.tsv', sep='\t', header=0)

# Read the ICA mixing matrix  (adjusting path as necessary)
ica_mix = pd.read_csv('../desc-ICA_mixing.tsv', sep='\t', header=0)

File ~/GitHub/multi-echo-data-analysis/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:873, in read_csv(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, skip_blank_lines, parse_dates, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, low_memory, memory_map, float_precision, storage_options, dtype_backend)
kwds_defaults = _refine_defaults_read(
   dialect,
   delimiter,
   (...)    869     dtype_backend=dtype_backend,
)
kwds.update(kwds_defaults)
--> 873 return _read(filepath_or_buffer, kwds)

File ~/GitHub/multi-echo-data-analysis/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:300, in _read(filepath_or_buffer, kwds)
_validate_names(kwds.get("names", None))
# Create the parser.
--> 300 parser = TextFileReader(filepath_or_buffer, **kwds)
if chunksize or iterator:
   return parser

File ~/GitHub/multi-echo-data-analysis/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1645, in TextFileReader.__init__(self, f, engine, **kwds)
   self.options["has_index_names"] = kwds["has_index_names"]
self.handles: IOHandles | None = None
-> 1645 self._engine = self._make_engine(f, self.engine)

File ~/GitHub/multi-echo-data-analysis/.venv/lib/python3.12/site-packages/pandas/io/parsers/readers.py:1904, in TextFileReader._make_engine(self, f, engine)
   if "b" not in mode:
       mode += "b"
-> 1904 self.handles = get_handle(
   f,
   mode,
   encoding=self.options.get("encoding", None),
   compression=self.options.get("compression", None),
   memory_map=self.options.get("memory_map", False),
   is_text=is_text,
   errors=self.options.get("encoding_errors", "strict"),
   storage_options=self.options.get("storage_options", None),
)
assert self.handles is not None
f = self.handles.handle

File ~/GitHub/multi-echo-data-analysis/.venv/lib/python3.12/site-packages/pandas/io/common.py:926, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
elif isinstance(handle, str):
   # Check whether the filename is to be opened in binary mode.
   # Binary mode does not support 'encoding' and 'newline'.
   if ioargs.encoding and "b" not in ioargs.mode:
       # Encoding
--> 926         handle = open(
           handle,
           ioargs.mode,
           encoding=ioargs.encoding,
           errors=errors,
           newline="",
       )
   else:
       # Binary mode
       handle = open(handle, ioargs.mode)

FileNotFoundError: [Errno 2] No such file or directory: '../manual_classification.tsv'

need [TO POTENTIALLY EXTRACT NOISE AND NON-NOISE IC TIMESERIES AND] a brain mask and, if you want,

After installing phys2cvr, potentially with extra dependencies, we can import the main workflow of phys2cvr and call its help to see all available parameters (many).

from phys2cvr.workflows import phys2cvr as p2c

help(p2c)

You may notice many and one potential parameters. phys2cvr is meant to be very versatile, but fear not! Only few of them are needed in different situations. For more info, check the workflow API page

As default, phys2cvr takes anything given to it and creates lagged regressors of the interpolated end-tidal of a given trace (that is, the \(P_{ET}CO_2hrf\)) after aligning it with a reference signal and convolving it with a canonical HRF. However, you can run a lag regression internally. Generally speaking, you may need these parameters - although only fname_func is mandatory:

  fname_func='../optcom.nii.gz',      # The optimally combined data
  fname_co2='../co2.phys',            # The CO2 trace with peaks
  fname_mask='../brain_mask.nii.gz',  # A mask to limit the analysis to
  fname_roi='../roi_mask.nii.gz',     # A ROI to use as signal reference
  outdir='../outdir',                 # Output folder
  run_regression=True,                # Run the analysis!
  lag_max=9,                          # Maximum lag to consider
  lag_step=0.3,                       # Lag step to consider
  l_degree=2,                         # Legendre polynomials degrees

You can also run a simple regression instead of a lagged one - in that case, remove lag_max and lag_step.

In terms of regressor for your CVR analysis, many options are supported:

\(CO_2\) traces (and their peak indices)
pre-made \(P_{ET}CO_2\) traces
task designs (timeseries)
BOLD signals averages from a ROI mask.

\(CO_2\) traces (and their peak indices)#

You can get \(CO_2\) traces (and their peak indices) with whatever program you want - simply run a very light lowpass filter (and maybe a downsampling), followed by peak detection.

Here’s an example with Physiopy’s peakdet.

import os

import numpy as np
import peakdet as pk

co2_path = os.path.abspath('../co2data')

raise ValueError("SKIP")
data = np.genfromtxt(co2_path)

# Set the desired "channel", i.e. column, of the file containing CO2 data
channel=0
# Set the sampling frequency of the data
samp_freq=10000

co2 = pk.Physio(data[:, channel], fs=samp_freq)

# Interpolation (downsampling), then filtering.
co2 = pk.operations.interpolate_physio(co2, 40.0, kind='linear')
co2 = pk.operations.filter_physio(co2, 2, method='lowpass', order=7)

# Automatic peak detection, then manual editing
# If you see that not enough/too many peaks get selected,
# you might tune the threshold and distance.
co2 = pk.operations.peakfind_physio(co2, thresh=0.5, dist=120)
co2 = pk.operations.edit_physio(co2)

outfile = '../co2.phys'
path = pk.save_physio(outfile, co2)

With this last output, we can run phys2cvr!

p2c(
  fname_func='../optcom.nii.gz',      # The optimally combined data
  fname_co2='../co2.phys',            # The CO2 trace with peaks
  fname_mask='../brain_mask.nii.gz',  # A mask to limit the analysis to
  fname_roi='../roi_mask.nii.gz',     # A ROI to use as signal reference
  outdir='../outdir',                 # Output folder
  run_regression=True,                # Run the analysis!
  lag_max=9,                          # Maximum lag to consider
  lag_step=0.3,                       # Lag step to consider
  l_degree=2,                         # Legendre polynomials degrees
  denoise_matrix_file=['../motion_params', '../other_noise'],  # All the noise you want to add
  orthogonalised_matrix_file=['../rejected_ic_timeseries.csv'],
  extra_matrix_file=['../accepted_ic_timeseries.csv'],
  scale_factor=None,                  # If CO2 trace is not in mmHg, you can scale the final map accordingly
  )

If you didn’t use peakdet for your \(CO_2\) processing and peak detection, add the parameter fname_pidx='../co2_peak_indices to the list.

phys2cvr will compute the \(P_{ET}CO_2hrf\) trace by linearly interpolating the end-tidal peaks of your \(CO_2\) trace and convolving the former with a canonical HRF. You can select another response function (as well as skipping the interpolation) using the parameter response_function.

Pre-made \(P_{ET}CO_2\) traces#

If you have a pre-made \(P_{ET}CO_2\) trace, simply skip the steps to make one.

p2c(
  fname_func='../optcom.nii.gz',      # The optimally combined data
  fname_co2='../petco2_trace',        # The pre-made PetCO2 trace
  fname_mask='../brain_mask.nii.gz',  # A mask to limit the analysis to
  fname_roi='../roi_mask.nii.gz',     # A ROI to use as signal reference
  outdir='../outdir',                 # Output folder
  run_regression=True,                # Run the analysis!
  lag_max=9,                          # Maximum lag to consider
  lag_step=0.3,                       # Lag step to consider
  l_degree=2,                         # Legendre polynomials degrees
  denoise_matrix_file=['../motion_params', '../other_noise'],  # All the noise you want to add
  orthogonalised_matrix_file=['../rejected_ic_timeseries.csv'],
  extra_matrix_file=['../accepted_ic_timeseries.csv'],
  scale_factor=None,                  # If CO2 trace is not in mmHg, you can scale the final map accordingly
  comp_endtidal=False,                # DO NOT compute PetCO2 internally
  response_function=None,             # DO NOT interpolate the PetCO2 trace.
  )

Task designs (timeseries)#

Simply use your task design as \(CO_2\) trace. You can interpolate it with an HRF or a respiratory response function internally, or do this beforehand.

p2c(
  fname_func='../optcom.nii.gz',      # The optimally combined data
  fname_co2='../task_design',        # The pre-made PetCO2 trace
  fname_mask='../brain_mask.nii.gz',  # A mask to limit the analysis to
  fname_roi='../roi_mask.nii.gz',     # A ROI to use as signal reference
  outdir='../outdir',                 # Output folder
  run_regression=True,                # Run the analysis!
  lag_max=9,                          # Maximum lag to consider
  lag_step=0.3,                       # Lag step to consider
  l_degree=2,                         # Legendre polynomials degrees
  denoise_matrix_file=['../motion_params', '../other_noise'],  # All the noise you want to add
  orthogonalised_matrix_file=['../rejected_ic_timeseries.csv'],
  extra_matrix_file=['../accepted_ic_timeseries.csv'],
  scale_factor=None,                  # If CO2 trace is not in mmHg, you can scale the final map accordingly
  comp_endtidal=False,                # DO NOT compute end-tidal trace
  response_function='rrf',            # Interpolate with RRF
  )

BOLD signals averages from a ROI mask.#

This is the simplest way to run phys2cvr, using the average BOLD signal from a desired ROI.

p2c(
  fname_func='../optcom.nii.gz',      # The optimally combined data
  fname_mask='../brain_mask.nii.gz',  # A mask to limit the analysis to
  fname_roi='../roi_mask.nii.gz',     # A ROI to use as signal reference
  outdir='../outdir',                 # Output folder
  run_regression=True,                # Run the analysis!
  lag_max=9,                          # Maximum lag to consider
  lag_step=0.3,                       # Lag step to consider
  l_degree=2,                         # Legendre polynomials degrees
  denoise_matrix_file=['../motion_params', '../other_noise'],  # All the noise you want to add
  orthogonalised_matrix_file=['../rejected_ic_timeseries.csv'],
  extra_matrix_file=['../accepted_ic_timeseries.csv'],
  comp_endtidal=False,                # DO NOT compute end-tidal trace
  response_function=None,             # DO NOT interpolate the trace
  apply_filter=True,                  # Filter the average BOLD signal (optional)
  highcut=0.04,                       # Lowpass frequency
  lowcut=0.02,                        # Highpass frequency
  butter_order=9,                     # Butterworth filter order
  )

Final remarks#

Remember to cite the right papers and resources depending on the desired pipeline!

For more info about phys2cvr, as well as knowing what to cite in which case, check out its documentation.

Not a python person? Prefer to run bash CLI commands? Not a problem. phys2cvr is in fact primarily designed to be used via bash CLI

phys2cvr --help

For more info about ME CVR, check out Moia et al. [2021], Cohen and Wang [2019], and .

If you want to contribute to phys2cvr, have requests, comments, or questions, feel free to open an issue!