Adding Synthetic Hemodynamic Reponses to Data
This example notebook illustrates the functionality in cedalion.sim.synthetic_hrf to create simulated datasets with added activations.
[1]:
# set this flag to True to enable interactive 3D plots
INTERACTIVE_PLOTS = False
[2]:
import os
import matplotlib.pyplot as plt
import numpy as np
import pyvista as pv
import xarray as xr
import cedalion
import cedalion.dataclasses as cdc
import cedalion.datasets
import cedalion.geometry.landmarks as cd_landmarks
import cedalion.imagereco.forward_model as fw
import cedalion.models.glm as glm
import cedalion.nirs
import cedalion.plots
import cedalion.sigproc.quality as quality
import cedalion.sim.synthetic_hrf as synhrf
import cedalion.xrutils as xrutils
from cedalion import units
from cedalion.imagereco.solver import pseudo_inverse_stacked
xr.set_options(display_expand_data=False)
if INTERACTIVE_PLOTS:
pv.set_jupyter_backend('server')
else:
pv.set_jupyter_backend('static')
Loading and preprocessing the dataset
This notebook uses a high-density, whole head resting state dataset recorded with a NinjaNIRS 22.
[3]:
rec = cedalion.datasets.get_nn22_resting_state()
geo3d = rec.geo3d
meas_list = rec._measurement_lists["amp"]
amp = rec["amp"]
amp = amp.pint.dequantify().pint.quantify("V")
display(amp)
Downloading file 'nn22_resting_state.zip' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/nn22_resting_state.zip' to '/home/runner/.cache/cedalion'.
Unzipping contents of '/home/runner/.cache/cedalion/nn22_resting_state.zip' to '/home/runner/.cache/cedalion/nn22_resting_state.zip.unzip'
<xarray.DataArray (channel: 567, wavelength: 2, time: 3309)> Size: 30MB
[V] 0.0238 0.02385 0.02375 0.02361 0.02364 ... 0.2887 0.2892 0.2897 0.291 0.2923
Coordinates:
* time (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0
samples (time) int64 26kB 0 1 2 3 4 5 ... 3303 3304 3305 3306 3307 3308
* channel (channel) object 5kB 'S10D87' 'S10D94' ... 'S47D29' 'S5D137'
source (channel) object 5kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (channel) object 5kB 'D87' 'D94' 'D89' ... 'D129' 'D29' 'D137'
* wavelength (wavelength) float64 16B 760.0 850.0
Attributes:
data_type_group: unprocessed raw[4]:
cedalion.plots.plot_montage3D(rec["amp"], geo3d)
[5]:
# Select channels which have at least a signal-to-noise ratio of 10
snr_thresh = 10 # the SNR (std/mean) of a channel.
snr, snr_mask = quality.snr(rec["amp"], snr_thresh)
amp_selected, masked_channels = xrutils.apply_mask(
rec["amp"], snr_mask, "drop", "channel"
)
print(f"Removed {len(masked_channels)} channels because of low SNR.")
Removed 22 channels because of low SNR.
[6]:
# Calculate optical density
od = cedalion.nirs.int2od(amp_selected)
Construct headmodel
We load the the Colin27 headmodel, since we need the geometry for image reconstruction.
[7]:
SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_colin27_segmentation()
[8]:
head_ijk = fw.TwoSurfaceHeadModel.from_surfaces(
segmentation_dir=SEG_DATADIR,
mask_files = mask_files,
brain_surface_file= os.path.join(SEG_DATADIR, "mask_brain.obj"),
scalp_surface_file= os.path.join(SEG_DATADIR, "mask_scalp.obj"),
landmarks_ras_file=landmarks_file,
brain_face_count=None,
scalp_face_count=None,
fill_holes=True, # needs to be true, otherwise landmark calculation fails
)
# transform coordinates to a RAS coordinate system
head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)
[9]:
display(head_ijk.brain)
display(head_ras.brain)
TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(15002, 3), faces.shape=(29988, 3))>, crs='ijk', units=<Unit('dimensionless')>)
TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(15002, 3), faces.shape=(29988, 3))>, crs='aligned', units=<Unit('millimeter')>)
[10]:
head_ras.landmarks
[10]:
<xarray.DataArray (label: 4, aligned: 3)> Size: 96B
[mm] -0.04992 79.79 -36.14 1.756 -100.9 ... -16.12 -54.66 75.95 -13.04 -54.08
Coordinates:
* label (label) <U3 48B 'Nz' 'Iz' 'LPA' 'RPA'
type (label) object 32B PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: alignedhead.landmarks contains the 4 landmarks [‘Nz’ ‘Iz’ ‘LPA’ ‘RPA’]. Since we want to create synthetic HRFs on the brain surface at landmark positions, we need to build the remaining 10-10 landmarks
[11]:
lmbuilder = cd_landmarks.LandmarksBuilder1010(head_ras.scalp, head_ras.landmarks)
all_landmarks = lmbuilder.build()
head_ras.landmarks = all_landmarks
/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:242: UserWarning: WIP: distance calculation around ears
warnings.warn("WIP: distance calculation around ears")
[12]:
center_brain = np.mean(head_ras.brain.mesh.vertices, axis=0)
We want to build the synthetic HRFs at C3 and C4 (green dots in the image below)
[13]:
plt_pv = pv.Plotter()
cedalion.plots.plot_surface(plt_pv, head_ras.brain, color="#d3a6a1")
cedalion.plots.plot_surface(plt_pv, head_ras.scalp, opacity=0.1)
cedalion.plots.plot_labeled_points(
plt_pv, head_ras.landmarks.sel(label=["C3", "C4"]), show_labels=True
)
plt_pv.camera.position = (-400, 500,400)
plt_pv.show()
Build spatial activation pattern on brain surface for landmarks C3 and C4
The function build_spatial_activation is used to place a spatial activation pattern on the brain surface. The activation pattern is a Gaussian function of the geodesic distance to a seed vertex. Hence, the size of the activation is determined by the standard deviation of this Gaussian, specified by the parameter spatial_scale. The peak intensity in HbO is determined by the parameter intensity_scale. The intensity of HbR activation is specified relative to the HbO peak intensity. So
if the HbO pattern describes an increase in Hbo then providing a negative factor smaller than 1 yields a decrease in HbR with smaller amplitude. The seed vertex (integer) can be selected as the closest vertex to a given landmark:
[14]:
# obtain the closest vertices to C3 and C4
c3_seed = head_ras.brain.mesh.kdtree.query(head_ras.landmarks.sel(label='C3'))[1]
c4_seed = head_ras.brain.mesh.kdtree.query(head_ras.landmarks.sel(label='C4'))[1]
/home/runner/miniconda3/envs/cedalion/lib/python3.11/site-packages/xarray/core/variable.py:338: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
data = np.asarray(data)
[15]:
# create the spatial activation
spatial_act = synhrf.build_spatial_activation(
head_ras,
c3_seed,
spatial_scale=2 * cedalion.units.cm,
intensity_scale=1 * units.micromolar,
hbr_scale=-0.4,
)
/home/runner/miniconda3/envs/cedalion/lib/python3.11/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:603: SparseEfficiencyWarning: splu converted its input to CSC format
return splu(A).solve
The resulting DataArray contains an activation value for each vertex and chromophore on the brain surface.
[16]:
display(spatial_act)
<xarray.DataArray (vertex: 15002, chromo: 2)> Size: 240kB [M] 1.168e-19 -4.672e-20 1.648e-20 -6.593e-21 3.798e-21 ... 0.0 -0.0 0.0 -0.0 Coordinates: * chromo (chromo) <U3 24B 'HbO' 'HbR' Dimensions without coordinates: vertex
[17]:
f,ax = plt.subplots(1,2,figsize=(10,5))
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_act.sel(chromo="HbO").pint.to("uM"),
head_ras.brain,
ax[0],
camera_pos="C3",
cmap="RdBu_r",
vmin=-1,
vmax=+1,
cb_label=r"$\Delta$ HbO / µM",
title=None,
)
ax[0].set_title("HbO")
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_act.sel(chromo="HbR").pint.to("uM"),
head_ras.brain,
ax[1],
camera_pos="C3",
cmap="RdBu_r",
vmin=-1,
vmax=+1,
cb_label=r"$\Delta$ HbR / µM",
title=None,
)
ax[1].set_title("HbR");
/home/runner/miniconda3/envs/cedalion/lib/python3.11/site-packages/xarray/core/variable.py:338: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
data = np.asarray(data)
The following plot illustrates the effects of the spatial_scale and intensity_scale parameters:
[18]:
f, ax = plt.subplots(2, 3, figsize=(9, 6))
for i, spatial_scale in enumerate([0.5 * units.cm, 2 * units.cm, 3 * units.cm]):
spatial_act = synhrf.build_spatial_activation(
head_ras,
c3_seed,
spatial_scale=spatial_scale,
intensity_scale=1 * units.micromolar,
hbr_scale=-0.4,
)
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_act.sel(chromo="HbO").pint.to("uM"),
head_ras.brain,
ax[0, i],
camera_pos="C3",
cmap="RdBu_r",
vmin=-1,
vmax=+1,
cb_label=r"$\Delta$ HbO / µM",
title=None,
)
ax[0, i].set_title(f"spatial_scale: {spatial_scale.magnitude} cm")
for i, intensity_scale in enumerate(
[
0.5 * units.micromolar,
1.0 * units.micromolar,
2.0 * units.micromolar,
]
):
spatial_act = synhrf.build_spatial_activation(
head_ras,
c3_seed,
spatial_scale=2 * units.cm,
intensity_scale=intensity_scale,
hbr_scale=-0.4,
)
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_act.sel(chromo="HbO").pint.to("uM"),
head_ras.brain,
ax[1, i],
camera_pos="C3",
cmap="RdBu_r",
vmin=-2,
vmax=+2,
cb_label=r"$\Delta$ HbO / µM",
title=None,
)
ax[1, i].set_title(f"intensity_scale: {intensity_scale.magnitude} µM")
f.tight_layout()
For this example notebook two activations are placed below C3 and C4:
[19]:
spatial_act_c3 = synhrf.build_spatial_activation(
head_ras,
c3_seed,
spatial_scale=2 * cedalion.units.cm,
intensity_scale=1 * units.micromolar,
hbr_scale=-0.4,
)
spatial_act_c4 = synhrf.build_spatial_activation(
head_ras,
c4_seed,
spatial_scale=2 * cedalion.units.cm,
intensity_scale=1 * units.micromolar,
hbr_scale=-0.4,
)
We concatenate the two images for C3 and C4 along dimension trial_type to get a single DataArray with the spatial information for both landmarks.
[20]:
# concatenate the two spatial activations along a new dimension
spatial_imgs = xr.concat(
[spatial_act_c3, spatial_act_c4], dim="trial_type"
).assign_coords(trial_type=["Stim C3", "Stim C4"])
spatial_imgs
[20]:
<xarray.DataArray (trial_type: 2, vertex: 15002, chromo: 2)> Size: 480kB [M] 1.168e-19 -4.672e-20 1.648e-20 -6.593e-21 ... -2.275e-13 1.01e-19 -4.038e-20 Coordinates: * chromo (chromo) <U3 24B 'HbO' 'HbR' * trial_type (trial_type) <U7 56B 'Stim C3' 'Stim C4' Dimensions without coordinates: vertex
Plots of spatial patterns
Using the helper function cedalion.plots.brain_plot, the created activations on the brain surface below C3 and C4 are plotted:
[21]:
f, ax = plt.subplots(1,2, figsize=(10,5))
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_imgs.sel(trial_type="Stim C3", chromo="HbO").pint.to("uM"),
head_ras.brain,
ax[0],
camera_pos="C3",
cmap="RdBu_r",
vmin=-1,
vmax=+1,
cb_label=r"$\Delta$ HbO / µM",
title="C3 Activation",
)
cedalion.plots.brain_plot(
od,
head_ras.landmarks,
spatial_imgs.sel(trial_type="Stim C4", chromo="HbO").pint.to("uM"),
head_ras.brain,
ax[1],
camera_pos="C4",
cmap="RdBu_r",
vmin=-1,
vmax=+1,
cb_label=r"$\Delta$ HbO / µM",
title="C4 Activation",
)
f.tight_layout()
Image Reconstruction
We load the precomputed Adot matrix to be able to map from image to channel space. (For details see image_reconstruction example notebook).
[22]:
Adot = cedalion.datasets.get_ninjanirs_colin27_precomputed_sensitivity()
[23]:
# we only consider brain vertices, not scalp
Adot_brain = Adot[:, (Adot.is_brain).values,:]
# drop the pruned channels
Adot_brain = Adot_brain.sel(channel=od.channel)
Adot_brain
[23]:
<xarray.DataArray (channel: 545, vertex: 15002, wavelength: 2)> Size: 131MB
[16352180 values with dtype=float64]
Coordinates:
is_brain (vertex) bool 15kB True True True True ... True True True True
* channel (channel) object 4kB 'S10D87' 'S10D94' ... 'S47D29' 'S5D137'
* wavelength (wavelength) float64 16B 760.0 850.0
source (channel) object 4kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (channel) object 4kB 'D87' 'D94' 'D89' ... 'D129' 'D29' 'D137'
Dimensions without coordinates: vertexThe forward model and image reconstruction translate between timeseries of different wavelengths in channel space and time series of different chromophores in image space. To this end the image reconstruction operates on stacked arrays in which the dimensions ‘channel’ and ‘wavelength’ are stacked to form a new dimension ‘flat_channel’. Likewise the dimensions ‘vertex’ and ‘chromo’ are stacked as ‘flat_vertex’.
[24]:
Adot_stacked = fw.ForwardModel.compute_stacked_sensitivity(Adot_brain)
Adot_stacked = Adot_stacked.pint.quantify()
Adot_stacked
[24]:
<xarray.DataArray (flat_channel: 1090, flat_vertex: 30004)> Size: 262MB
[1/M] 5.859e-16 4.574e-15 3.6e-15 2.592e-15 ... 1.471e-11 0.0108 1.169e-10
Coordinates:
is_brain (flat_vertex) bool 30kB True True True True ... True True True
chromo (flat_vertex) <U3 360kB 'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR' 'HbR'
vertex (flat_vertex) int64 240kB 0 1 2 3 4 ... 14998 14999 15000 15001
wavelength (flat_channel) float64 9kB 760.0 760.0 760.0 ... 850.0 850.0
channel (flat_channel) object 9kB 'S10D87' 'S10D94' ... 'S5D137'
source (flat_channel) object 9kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (flat_channel) object 9kB 'D87' 'D94' 'D89' ... 'D29' 'D137'
Dimensions without coordinates: flat_channel, flat_vertexinvert the sensitivity matrix:
[25]:
Adot_inverted = pseudo_inverse_stacked(Adot_stacked)
Adot_inverted = Adot_inverted.pint.quantify()
Adot_inverted
[25]:
<xarray.DataArray (flat_vertex: 30004, flat_channel: 1090)> Size: 262MB
[M] 2.581e-13 -7.848e-13 -5.654e-12 ... -2.363e-14 2.218e-14 -4.732e-15
Coordinates:
is_brain (flat_vertex) bool 30kB True True True True ... True True True
chromo (flat_vertex) <U3 360kB 'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR' 'HbR'
vertex (flat_vertex) int64 240kB 0 1 2 3 4 ... 14998 14999 15000 15001
wavelength (flat_channel) float64 9kB 760.0 760.0 760.0 ... 850.0 850.0
channel (flat_channel) object 9kB 'S10D87' 'S10D94' ... 'S5D137'
source (flat_channel) object 9kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (flat_channel) object 9kB 'D87' 'D94' 'D89' ... 'D29' 'D137'
Dimensions without coordinates: flat_vertex, flat_channelTo multiply the spatial image with the sensitivity matrix the spatial image’s vertex and chromo dimensions must be stacked, too.
[26]:
spatial_imgs_stacked = fw.stack_flat_vertex(spatial_imgs)
display(spatial_imgs_stacked)
<xarray.DataArray (trial_type: 2, flat_vertex: 30004)> Size: 480kB [M] 1.168e-19 1.648e-20 3.798e-21 1.274e-19 ... -3.851e-12 -2.275e-13 -4.038e-20 Coordinates: * trial_type (trial_type) <U7 56B 'Stim C3' 'Stim C4' * flat_vertex (flat_vertex) object 240kB MultiIndex * chromo (flat_vertex) <U3 360kB 'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR' 'HbR' * vertex (flat_vertex) int64 240kB 0 1 2 3 4 ... 14998 14999 15000 15001
We can now map our spatial patterns to channel space:
[27]:
spatial_chan_stacked = Adot_stacked @ spatial_imgs_stacked
spatial_chan_stacked
[27]:
<xarray.DataArray (flat_channel: 1090, trial_type: 2)> Size: 17kB
[] -1.963e-16 -6.505e-07 -5.609e-16 -6.902e-08 ... 1.986e-07 3.311e-15 7.887e-10
Coordinates:
wavelength (flat_channel) float64 9kB 760.0 760.0 760.0 ... 850.0 850.0
channel (flat_channel) object 9kB 'S10D87' 'S10D94' ... 'S5D137'
source (flat_channel) object 9kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (flat_channel) object 9kB 'D87' 'D94' 'D89' ... 'D29' 'D137'
* trial_type (trial_type) <U7 56B 'Stim C3' 'Stim C4'
Dimensions without coordinates: flat_channel[28]:
spatial_chan = fw.unstack_flat_channel(spatial_chan_stacked)
Show the spatial activation in channel space with a scalp plot.
[29]:
fig, ax = plt.subplots(1, 2)
# adjust plot size
fig.set_size_inches(12, 6)
cedalion.plots.scalp_plot(
od,
rec.geo3d,
spatial_chan.sel(trial_type="Stim C3", wavelength=850),
ax[0],
cmap="YlOrRd",
title="850nm, activation under C3",
vmin=spatial_chan.values.min(),
vmax=spatial_chan.values.max(),
cb_label="max peak amplitude",
)
cedalion.plots.scalp_plot(
od,
rec.geo3d,
spatial_chan.sel(trial_type="Stim C4", wavelength=850),
ax[1],
cmap="YlOrRd",
title="850nm, activation under C4",
vmin=spatial_chan.values.min(),
vmax=spatial_chan.values.max(),
cb_label="Max peak amplitude",
)
plt.show()
Get top 5 channels for each trial type where synthetic activation is highest
[30]:
roi_chans_c3 = spatial_chan.channel[
spatial_chan.sel(trial_type="Stim C3").max("wavelength").argsort()[-5:].values
].values
roi_chans_c4 = spatial_chan.channel[
spatial_chan.sel(trial_type="Stim C4").max("wavelength").argsort()[-5:].values
].values
Concentration Scale
The activations were simulataed in image space with a peak concentration change of 1 µM in one vertex. The change in optical density in one channel reflects concentration changes in the ensemble of vertices that this channel is sensitive to.
When applying the Beer-Lambert-transformation in channel space, a change in concentration is calculated for each channel. However, the scales of these concentration changes and the concentration changes in single vertices are not the same.
Here, a correction factor is calculated to scale the activation in channel space to 1uM.
[31]:
dpf = xr.DataArray(
[6, 6],
dims="wavelength",
coords={"wavelength": rec["amp"].wavelength},
)
[32]:
# add time axis with one time point so we can convert to conc
spatial_chan_w_time = spatial_chan.expand_dims("time")
spatial_chan_w_time = spatial_chan_w_time.assign_coords(time=[0])
spatial_chan_w_time.time.attrs["units"] = "second"
display(spatial_chan_w_time)
spatial_chan_conc = cedalion.nirs.od2conc(
spatial_chan_w_time, geo3d, dpf, spectrum="prahl"
)
<xarray.DataArray (time: 1, trial_type: 2, wavelength: 2, channel: 545)> Size: 17kB
[] -1.738e-14 -6.175e-13 -3.619e-13 -2.301e-15 ... 3.903e-06 3.561e-08 2.022e-08
Coordinates:
* wavelength (wavelength) float64 16B 760.0 850.0
* channel (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
* trial_type (trial_type) <U7 56B 'Stim C3' 'Stim C4'
source (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
detector (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
* time (time) int64 8B 0[33]:
# rescale so that synthetic hrfs add 1 micromolar at peak.
rescale_factor = (1* units.micromolar / spatial_chan_conc.max())
display(rescale_factor)
spatial_chan *= rescale_factor
<xarray.DataArray 'concentration' ()> Size: 8B [] 17.98
HRFs in channel space
So far the notebook focused on the spatial extent of the activation. To build the temporal HRF model we use the same functionality that generates hrf regressors for the GLM.
First we select a basis function, which defines the temporal shape of the HRF.
[34]:
basis_fct = glm.Gamma(tau=0 * units.s, sigma=3 * units.s, T=3 * units.s)
[35]:
od.time
[35]:
<xarray.DataArray 'time' (time: 3309)> Size: 26kB
0.0 0.1112 0.2225 0.3337 0.445 0.5562 ... 367.5 367.6 367.7 367.8 367.9 368.0
Coordinates:
* time (time) float64 26kB 0.0 0.1112 0.2225 0.3337 ... 367.8 367.9 368.0
samples (time) int64 26kB 0 1 2 3 4 5 6 ... 3303 3304 3305 3306 3307 3308
Attributes:
units: secondA Stim DataFrame, which contains the onset, duration and amplitude of the synthetic HRFs, is created.
[36]:
stim_df = synhrf.build_stim_df(
max_time=od.time.values[-1] * units.seconds,
trial_types=["Stim C3", "Stim C4"],
min_interval=10 * units.seconds,
max_interval=20 * units.seconds,
min_stim_dur = 10 * units.seconds,
max_stim_dur = 10 * units.seconds,
min_stim_value = 1.0,
max_stim_value = 1.0,
order="random",
)
[37]:
stim_df.head()
[37]:
| onset | duration | value | trial_type | |
|---|---|---|---|---|
| 0 | 17.35 | 10.0 | 1.0 | Stim C3 |
| 1 | 44.71 | 10.0 | 1.0 | Stim C4 |
| 2 | 71.19 | 10.0 | 1.0 | Stim C4 |
| 3 | 91.41 | 10.0 | 1.0 | Stim C4 |
| 4 | 117.11 | 10.0 | 1.0 | Stim C4 |
We can now use our stim dataframe, basis function, and spatial information to create the synthetic HRF timeseries
[38]:
syn_ts = synhrf.build_synthetic_hrf_timeseries(od, stim_df, basis_fct, spatial_chan)
We get a synthetic HRF timeseries for each channel, trial_type and chromo / wavelength
[39]:
syn_ts
[39]:
<xarray.DataArray (trial_type: 2, wavelength: 2, channel: 545, time: 3309)> Size: 58MB
[] -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Coordinates:
* wavelength (wavelength) float64 16B 760.0 850.0
* channel (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
* trial_type (trial_type) <U7 56B 'Stim C3' 'Stim C4'
source (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
detector (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
* time (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0Sum the synthetic timeseries over trial_type dimension, so it has the same shape as the resting state data
[40]:
syn_ts_sum = syn_ts.sum(dim='trial_type')
syn_ts_sum
[40]:
<xarray.DataArray (wavelength: 2, channel: 545, time: 3309)> Size: 29MB
[] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
Coordinates:
* wavelength (wavelength) float64 16B 760.0 850.0
* channel (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
source (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
detector (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
* time (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0Adding HRFs to measured data
Here, the simulated activations are combined with physiological noise by adding the synthetic HRFs to the resting state dataset:
[41]:
od_w_hrf = od + syn_ts_sum
Recover the HRFs again
In the following, the added activations should be extracted from the simulated dataset again. To this end, the data is frequency filtered and block averages are calculated.
[42]:
od_w_hrf_filtered = od_w_hrf.cd.freq_filter(fmin=0.02, fmax=0.5, butter_order=4)
[43]:
od_w_hrf_filtered
[43]:
<xarray.DataArray (channel: 545, wavelength: 2, time: 3309)> Size: 29MB
[] 0.002072 0.002634 0.00317 0.003662 ... 0.006999 0.006033 0.004973 0.00387
Coordinates:
* time (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0
* channel (channel) object 4kB 'S10D87' 'S10D94' ... 'S47D29' 'S5D137'
* wavelength (wavelength) float64 16B 760.0 850.0
samples (time) int64 26kB 0 1 2 3 4 5 ... 3303 3304 3305 3306 3307 3308
source (channel) object 4kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
detector (channel) object 4kB 'D87' 'D94' 'D89' ... 'D129' 'D29' 'D137'[44]:
epochs = od_w_hrf_filtered.cd.to_epochs(
stim_df, # stimulus dataframe
["Stim C3", "Stim C4"], # select events
before=5 * units.seconds, # seconds before stimulus
after=20 * units.seconds, # seconds after stimulus
)
# calculate baseline
baseline = epochs.sel(reltime=(epochs.reltime < 0)).mean("reltime")
# subtract baseline
epochs_blcorrected = epochs - baseline
# group trials by trial_type. For each group individually average the epoch dimension
blockaverage = epochs_blcorrected.groupby("trial_type").mean("epoch")
n_roi = roi_chans_c3.size
# show results
f, ax = plt.subplots(2, n_roi, figsize=(16, 8))
ax = ax.flatten()
for i_ch, ch in enumerate(roi_chans_c3):
for ls, trial_type in zip(["-", "--"], blockaverage.trial_type):
ax[i_ch].plot(
blockaverage.reltime,
blockaverage.sel(wavelength=760, trial_type=trial_type, channel=ch),
"r",
lw=2,
ls=ls,
)
ax[i_ch].plot(
blockaverage.reltime,
blockaverage.sel(wavelength=850, trial_type=trial_type, channel=ch),
"b",
lw=2,
ls=ls,
)
ax[i_ch].grid(1)
ax[i_ch].set_title(ch)
ax[i_ch].set_ylim(-0.05, 0.05)
for i_ch, ch in enumerate(roi_chans_c4):
for ls, trial_type in zip(["-", "--"], blockaverage.trial_type):
ax[i_ch + n_roi].plot(
blockaverage.reltime,
blockaverage.sel(wavelength=760, trial_type=trial_type, channel=ch),
"r",
lw=2,
ls=ls,
)
ax[i_ch + n_roi].plot(
blockaverage.reltime,
blockaverage.sel(wavelength=850, trial_type=trial_type, channel=ch),
"b",
lw=2,
ls=ls,
)
ax[i_ch + n_roi].grid(1)
ax[i_ch + n_roi].set_title(ch)
ax[i_ch + n_roi].set_ylim(-0.05, 0.05)
plt.suptitle(
"Blockaverage for channels most sensitive to C3 (top) and C4 (bottom): 760nm: r | 850nm: b | C3: - | C4: --"
)
plt.tight_layout()
plt.show()
Map block average back to brain surface
We map our extracted block averages back to the brain surface to visualize the recovered HRFs activation for Stim C3. We can compare it to the synthetic HRF image we created earlier.
[45]:
blockaverage_img = Adot_inverted @ fw.stack_flat_channel(blockaverage)
[46]:
# build an xindex to use .sel along the chromo dimension
blockaverage_img = blockaverage_img.set_xindex("chromo")
[47]:
# plot HbO time trace of left and right brain hemisphere during FTapping/Right
for view in ["left_hemi", "right_hemi"]:
trial_type = "Stim C3"
gif_fname = "Ftapping-right" + "_HbO_" + view + ".gif"
hbo = (
blockaverage_img.sel(chromo="HbO", trial_type=trial_type).pint.dequantify()
/ 1e-5
) # FIXME unit handling
hbo_brain = hbo
ntimes = hbo.sizes["reltime"]
b = cdc.VTKSurface.from_trimeshsurface(head_ras.brain)
b = pv.wrap(b.mesh)
b["reco_hbo"] = hbo_brain[:, 0] - hbo_brain[:, 0]
p = pv.Plotter()
p.add_mesh(
b,
scalars="reco_hbo",
cmap="seismic", # 'gist_earth_r',
clim=(-2.5, 2.5),
scalar_bar_args={"title": "HbO / µM"},
smooth_shading=True,
)
tl = lambda tt: f"{trial_type} HbO rel. time: {tt:.3f} s"
time_label = p.add_text(tl(0))
cog = head_ras.brain.vertices.mean("label").values
if view == "left_hemi":
p.camera.position = cog + [-400, 0, 0]
else:
p.camera.position = cog + [400, 0, 0]
p.camera.focal_point = cog
p.camera.up = [0, 0, 1]
p.reset_camera()
p.open_gif(gif_fname)
for i in range(0, ntimes, 3):
b["reco_hbo"] = hbo_brain[:, i] - hbo_brain[:, 0]
time_label.set_text("upper_left", tl(hbo_brain.reltime[i]))
p.write_frame()
p.close()
[48]:
from IPython.display import Image
display(Image(data=open("Ftapping-right_HbO_left_hemi.gif",'rb').read(), format='png'))
display(Image(data=open("Ftapping-right_HbO_right_hemi.gif",'rb').read(), format='png'))
