Calculating heart rate for each slice

[ckelly5]

Hi Lars,

I'm interested to calculate and export heart rate (HR) with more measurements than in physio.ons_secs.hr (which is HR per image volume). You previously gave me the helpful suggestion to try to calculate HR for each slice by setting the onset_slice to an array. I've tried doing this and I just want to check the steps and outputs. I'm running the PhysIO toolbox using matlab only with siemens_vd_ppu3t_matlab_script.m. I changed the onset slice by changing one line in siemens_vd_ppu3t_matlab_script.m, as follows:

Original: physio.scan_timing.sqpar.onset_slice = 1;
My change: physio.scan_timing.sqpar.onset_slice = [1:physio.scan_timing.sqpar.Nslices];

Then I ran siemens_vd_ppu3t_matlab_script.m, and it finished running without any errors, and the output was one physio.mat file for each slice. I have 72 slices in total, so the output is named physio_slice001.mat through to physio_slice072.mat. So I think that means I have HR for each slice (i.e., I have ons_secs.hr in each physio_slice.mat).

I can also change the sliding window width in tapas_physio_hr.m- that doesn't change the total number of HR measurements I get (I still get HR for each slice), but it changes the values of the HR measurements (I guess it is maybe changing how many measurements we're averaging over).

My question is- does that all sound correct and like it's working like you would expect? Do the outputs sound like what you would expect?

If so, I also was hoping to check the ordering of the HR measurements. In each of those physio_slice001.mat through to physio_slice072.mat, ons_secs.hr is 242x1. I have 242 volumes, so I assume ons_secs.hr for physio_slice001.mat is the HR of the first slice of each volume, and then ons_secs.hr for physio_slice002.mat is the HR of the second slice of each volume, etc. Is that right? If so, I guess that is ordered differently to how it's acquired (i.e. we would have acquired slice 1-72 for volume 1, then slice 1-72 for volume 2, etc, rather than slice 1 for volume 1-242, then slice 2 for volume 1-242, etc). So I just need to keep this in mind when I'm using the HR measurements, and potentially do a few extra steps to reorder the HR measurements if needed (e.g. to line up with condition timings).

I would really appreciate your help.
Many thanks in advance,
Claire

[mrikasper]

Dear Claire,

Thank you for trying out this functionality and giving such a detailed report.

Your strategy is perfectly sound, and I apologize that the PhysIO interface is clunky for your use case.

You are right that the physio_slice001.mat...physio_slice072.mat all contain one value per volume, sampled at the respective onset slice time. And onset slice here is meant temporally (not spatially where you would have to consider interleaved slice order/multi-band etc.), so the data in the file labeled slice001 is indeed sampled at time 0 within each respective TR, slice_002 at TR/nSlices, slice_003 at 2*TR/nSlices etc.

So you can just concatenate the different ons_secs.hr and have to interleave them, as you describe:

If so, I guess that is ordered differently to how it's acquired (i.e. we would have acquired slice 1-72 for volume 1, then slice 1-72 for volume 2, etc, rather than slice 1 for volume 1-242, then slice 2 for volume 1-242, etc)

If ons_secs_slice001.hr and ons_secs_slice002.hr are both 242x1, something like that should do the interleaving:

hr = [];
for iSlice = 1:nSlices
    load(sprintf('physio_slice%03d.mat', iSlice));
    hr[:,iSlice] = ons_secs.hr;
end
hr = reshape(hr', [],1); % the ' is important to create the transpose and mix row/column vectors correctly

By the way, If you just use PhysIO for this computation, you could manipulate tapas_physio_main_create_regressors in line 335 (starting from a single onset slice in the script) to temporally behave like having multiple onset_slices

[convHRV, ons_secs.hr, verbose] = tapas_physio_create_hrv_regressors(...
onset_slice_tmp = sqpar.onset_slice;
sqpar.onset_slice = 1:sqpar.Nslices
                ons_secs, sqpar, model.hrv, verbose);
sqpar_onset_slice = onset_slice_tmp
%   [model.R, model.R_column_names] = append_regressors(model.R, model.R_column_names, ...
                convHRV, 'HR * CRF');

As you see, you would have to comment out appending the multiple regressors related to HRV (or correct the lengths of these vectors), but ons_secs.hr should be highly sampled and saved into physio.mat.

I hope that helps!

All the best,
Lars

[ckelly5]

Dear Lars,

Thanks for your reply. I'm glad to hear the steps and outputs I mentioned are sounding correct. Thanks also for the commands you provided for concatenating and interleaving the different ons_secs.hr for each slice- I've been able to try that first block of commands you provided, and they worked. After running those commands, I ended up with the hr variable, which looks correct from my checks (it is 17424x1 - which is 242 volumes by 72 slices, ordered slice 1-72 for volume 1, then slice 1-72 for volume 2, etc). My next step is to line up this hr variable with the conditions/blocks of our task-based fMRI sequence. I have the onset time and duration in seconds of the conditions, so I need to figure out how many hr measurements fit into each condition, so I just need to make sure I understand the order and timing of the hr measurements in that hr variable.

Thanks for clarifying that the ordering is temporally not spatially. This is a multiband sequence (MB factor 3) and I believe I can get the spatial order of the slices from information in the dicom header, and also from the *_Info.log that I provided as input to the PhysIO Toolbox. So from my understanding, the hr variable (17424x1) is ordered temporally slice 1-72 for volume 1, then slice 1-72 for volume 2, etc., but spatially this would be different- spatially this is slice 0, 24, 48, 11, 35, 59, etc, for volume 1, and then the same for subsequent volumes (as shown in the *_Info.log). Does that sound correct? I hope I'm understanding it correctly, but please let me know if anything I mentioned sounds unclear or incorrect.

If I correctly understand the ordering of the hr measurements, then I think I can figure out the timing of each hr measurement from the slice timings, which I can get from the json file from running dcm2niix on the task fMRI sequence (e.g., I know from the dcm2niix json file that the first group of 3 slices- temporal slices 1,2,3, spatial slices 0,24,48- start at 0 seconds, then the second group of 3 slices start at 0.0825 seconds, then the third group of 3 slices start at 0.1625 seconds, etc). And once I have the timing of each hr measurement, I can figure out how many hr measurements fit into each fMRI task condition/block. Does that sound like a reasonable way to go about doing this?

In the dcm2niix json file, I get the slice timings for 1 volume (where the first group of 3 slices start at 0 seconds and the last group of 3 slices start at 1.895 seconds). So I still have to figure out the start time of volume 2 and all subsequent volumes. Our TR is 2 seconds, so I thought it's probably ok to assume the second volume starts at 2 seconds, and the third volume starts at 4 seconds, etc. Does that seem ok? Because another option is that I could figure out the start time of the second volume by adding the start time of the last slice of the first volume (1.895) and the typical slice duration (0.0825) = 1.9775 seconds. These two options give fairly similar but slightly different results (in terms of lining up the hr measures with the conditions), so I just wanted to ask about it.

I hope my questions make sense, but please let me know if I can provide any more information. Any feedback if possible would be really helpful!
Kind regards,
Claire