PhysIO - issues processing Siemens Vb data

[bmdeen]

Hi,

I recently used PhysIO (release 3.3.0) to process some old fMRI/physiological data from a Siemens Vb system, using the "nominal" option for scan time syncing, and providing a DCM file with scan timing info (sharing a sample script below). I was able to get this to work, but ran into a few issues that required me to modify the PhysIO scripts, that I wanted to mention here:

1) In my data, the respiration belt and ECG signals had different sampling rates (50Hz and 400Hz). In tapas_physio_read_physlogfiles_siemens, equal sampling rates are assumed, and lines 268-278 matches time durations by zero padding the shorter vector, which doesn't work when sampling rates differ. I replaced this code with code that interpolate the lower SR signal to the higher sampling rate.

2) My data was properly aligned to scan onset, read from the first DCM image, by tapas_read_physlogfiles_siemens, which aligned the time point ons_secs.t=0 to scan onset. However, running tapas_physio_preprocess_phys_timeseries.m undid this syncing, as lines 112-116 set ons_secs.t(1)=0 for plotting purposes. I was able to fix this by removing these lines of code.

Code I'm running:

`% PhysIO parameter settings physio = tapas_physio_new(); physio.save_dir = outDir; physio.log_files.vendor = 'Siemens'; physio.log_files.cardiac = ecgPath; physio.log_files.respiration = respPath; % physio.log_files.sampling_interval = 0.0025; % In our data, this is .0025 for cardiac, .02 for resp physio.log_files.relative_start_acquisition = 0; physio.log_files.align_scan = 'first'; physio.log_files.scan_timing = dcmPath; physio.preproc.cardiac.modality = 'ECG'; physio.verbose.level = 2; physio.verbose.process_log = cell(0, 1); physio.verbose.fig_handles = zeros(0, 1); physio.verbose.fig_output_file = 'physio.fig'; physio.verbose.use_tabs = false; physio.scan_timing.sqpar.TR = 2; physio.scan_timing.sqpar.Ndummies = 0; physio.scan_timing.sqpar.Nscans = 300; physio.scan_timing.sqpar.Nslices = 54; physio.scan_timing.sqpar.onset_slice = 1:54; % Slices to define regressors for physio.scan_timing.sync.method = 'nominal'; physio.preproc.cardiac.filter.include = false; physio.preproc.cardiac.initial_cpulse_select.method = 'auto_matched'; physio.preproc.cardiac.initial_cpulse_select.max_heart_rate_bpm = 120; physio.preproc.cardiac.initial_cpulse_select.min = 0.4; physio.preproc.cardiac.posthoc_cpulse_select.method = 'off'; physio.model.orthogonalise = 'none'; physio.model.censor_unreliable_recording_intervals = false; physio.model.output_multiple_regressors = 'multiple_regressors.txt'; physio.model.output_physio = 'physio.mat'; physio.model.retroicor.include = true; physio.model.retroicor.order.c = 3; physio.model.retroicor.order.r = 4; physio.model.retroicor.order.cr = 1; physio.model.rvt.include = true; physio.model.rvt.delays = 0; physio.model.hrv.include = true; physio.model.hrv.delays = 0; physio.model.noise_rois.include = false; physio.model.movement.include = false; physio.model.other.include = false; physio.ons_secs.c_scaling = 1; physio.ons_secs.r_scaling = 1;

% Run physiological recording preprocessing and noise modeling physio = tapas_physio_main_create_regressors(physio);`

Cheers and thanks for writing this helpful toolbox,

Ben

[mrikasper]

Dear Ben,

thank you for your feedback and for diving in to the toolbox. Concerning the version you ran, do you mean R2020a-v7.3.0 (check via tapas_physio_version)?

1. In principle, PhysIO should be able to work with different sampling rates for cardiac and respiratory data (you would have to provide a vector as log_files.sampling_interval, i.e., =[0.0025 0.02] s in your case. Does this solve the issue?

   • There could still be some trouble with zero-filling if there are unequal lengths of recordings. Thus, it's good to let us know this is an important feature to you, so we'll try to get it more robust.

2. The removal of the onset time in ons_secs.t(1)=0 should indeed only affect the plotting, and its value should be stored in ons_secs.t_start. Did you see a difference in the created regressors or the "cutout" plot of the time series used for physiological modeling that led you to conclude it was not working properly?

Thank you for your help in resolving these issues!

All the best, Lars

[bmdeen]

Hi Lars,

Thanks for the response. Yes, I'm using v7.3.0.

Re. sampling rate, I just tried using a vector definition of sampling rate, but it doesn’t seem to solve the issue. ons_secs.t and ons_secs.dt are still defined based on whichever sampling rate is higher - ons_secs.t by tapas_physio_read_physlogfiles_siemens, and ons_secs.dt by tapas_physio_preprocess_phys_timeseries. To account for differing sampling rates, I think the scripts would either need to resample one of the signals, or have separate definitions of ons_secs.t and ons_secs.dt for cardiac/repiratory data.

Re. onset time, I did see a difference in the cutout plot and RETROICOR regressors before/after removing the code that sets ons_secs.t(1)=0. The issue emerges when running tapas_physio_create_scan_timing_nominal, which doesn’t take t_start into account.

Best, Ben

[mrikasper]

Dear Ben,

Would you be able to share some example data to reproduce this behavior? Unfortunately I don't have data with differing sampling rates at hand.

With respect to the timing - why did you choose scan timing nominal, if you had the dicom file available? Just for comparison or was there another problem?

Thank you for your help!

All the best, Lars

[bmdeen]

I'm happy to share sample physio data. Not sure if I can share DCMs - would have to check our policy on this, since it's animal data. If the physio data is sufficient, let me know how to share it.

Re. using the nominal option, I may have misunderstood the appropriate option in this case - I suppose I should be using 'scan_timing_log'?

Thanks, Ben

[mrikasper]

Dear Ben,

If the files are not too big (<20MB), you can just send me an e-mail (see the README of PhysIO). Otherwise, sending me an online storage link or even putting it as an upload here on GitHub is fine with me.

If you cannot share the DICOMs, could you maybe save the variables tStartDicom and tStopDicom as computed in l.100-109 of tapas_physio_read_physlogfiles_siemens?

    dicomHeader             = spm_dicom_headers(...
        fullfile(log_files.scan_timing));

    tStartScanDicom    = dicomHeader{1}.AcquisitionTime;

    % TODO: Include AcquisitionNumber? InstanceNumber?
    tStopScanDicom     = dicomHeader{1}.AcquisitionTime + ...
        dicomHeader{1}.RepetitionTime/1000;

For the scan timing computation, indeed scan_timing_log is the correct selection if the DICOM timing shall be used. The difference to nominal depends on how much longer your physological logfile records after your scan ends. See tapas_physio_create_scan_timing l. 78ff vs 91ff.

I hope that helps!

All the best, Lars