SIEMENS PMU logging synchronization

Dear Lars, thank you for providing a toolbox for noise control in fMRI data! In our actual data set, we face a problem with synchronization of the PMU-logging files .resp and .puls with the raw image files MR... Independent of which MR..- file is given in „log_scan_time“, the toolbox takes into account the same cutout region for physiological regressors. When „align_scan“ is first, the very beginnig of the physiological recordings is cut out for regressors, when „align_scan“ is last, the very end of the physiological recordings are taken into for regressor building. In our scanning session, we did five runs, localizer, T1, T2 and two functional T2. Thus, we would expect the cutout regions fort he T2-runs to be in the middle of the physiological recordings, as PMU-logging was gathered from before all scanning started and ended after the last T2*-run.

Do you have any hint how to synchronize the data properly?

Kind regards, Alex

Dear Alex,

thank you for providing a toolbox for noise control in fMRI data!

Happy you are trying it out!

In our actual data set, we face a problem with synchronization of the PMU-logging files .resp and .puls with the raw image files MR... Independent of which MR..- file is given in „log_scan_time“, the toolbox takes into account the same cutout region for physiological regressors. When „align_scan“ is first, the very beginning of the physiological recordings is cut out for regressors, when „align_scan“ is last, the very end of the physiological recordings are taken into for regressor building.

I would need some more details on how you set up PhysIO, maybe you can send me a batch? Specifically, it seems you are using the nominal scan synchronization, because the align_scan option is mostly designed for that purpose. If you have a DICOM-file of your first fMRI volume (of the run you try to model), you could enter that into log_files.scan_timing and set scan_timing.sync.method = 'scan_timing_log';.

In our scanning session, we did five runs, localizer, T1, T2 and two functional T2. Thus, we would expect the cutout regions fort he T2-runs to be in the middle of the physiological recordings, as PMU-logging was gathered from before all scanning started and ended after the last T2*-run.

So you have one very long logfile? The raw data should be plotted by PhysIO before it tries the synchronization. Does the plotted raw trace have the right length (I would assume more than an hour?).

I hope that helps and look forward to your response, Lars

Dear Lars, thank you for your quick response!

I would need some more details on how you set up PhysIO, maybe you can send me a batch?

%----------------------------------------------------------------------- % Job saved on 02-Apr-2019 20:56:05 by cfg_util (rev $Rev: 6460 $) % spm SPM - SPM12 (6470) % cfg_basicio BasicIO - Unknown %----------------------------------------------------------------------- matlabbatch{1}.spm.tools.physio.save_dir = {'D:\38\t2stern_2'}; matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens'; matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'E:\rm2\PMU\20181212_Blinova_Karina.puls'}; matlabbatch{1}.spm.tools.physio.log_files.respiration = {'E:\rm2\PMU\20181212_Blinova_Karina.resp'}; matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'D:\38\20181212_13h28_0010915217_dicom\0035958628\1.3.12.2.1107.5.2.43.67036.201812121403419088128313.0.0.0\MR.18273.939'}; matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025; matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0; matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first'; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 12; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = []; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 1.5; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 677; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 4; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = 0.125; matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = []; matlabbatch{1}.spm.tools.physio.scan_timing.sync.scan_timing_log = struct([]); matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU'; matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4; matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat'; matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]); matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt'; matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat'; matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none'; matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false; matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3; matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4; matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1; matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]); matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]); matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]); matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]); matlabbatch{1}.spm.tools.physio.model.other.no = struct([]); matlabbatch{1}.spm.tools.physio.verbose.level = 2; matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = ''; matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

If you have a DICOM-file of your first fMRI volume (of the run you try to model), you could enter that into log_files.scan_timing and set scan_timing.sync.method = 'scan_timing_log';

In the batch I have the line

matlabbatch{1}.spm.tools.physio.scan_timing.sync.scan_timing_log = struct([]);

which comes close to your proposal (scan_timing.sync.method = 'scan_timing_log';) but is not exactly the same.

So you have one very long logfile? The raw data should be plotted by PhysIO before it tries the synchronization. Does the plotted raw trace have the right length (I would assume more than an hour?).

The logfiles are plotted for the appropriate time (about 3400s).

raw

When I put

matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first';

in the batch above:

cutout_first

When I put matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';

cutout_last

The Dicom file in the batch

matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'D:\38\20181212_13h28_0010915217_dicom\0035958628\1.3.12.2.1107.5.2.43.67036.201812121403419088128313.0.0.0\MR.18273.939'};

is the first scan of our second T2*-run and should appear in the last third of the Physio-logfile.

Hope that helps finding the culprit, Alex

Dear Alex,

the batch looks fine, and you have indeed set the sync.method correctly. The discrepancy between my suggested and you actual command is that I was referring to the native Matlab physio-structure defined in tapas_physio_new(), and you were using the Batch Editor matlabbatch version of that structure. They don't have the same fields, but when I check the matlabbatch with spm_jobman('interactive', matlabbatch), everything looks good.

So I am at a loss what might have happened, because the screenshots you kindly provided look as if the DICOM-timing is completely ignored. Also, the time axes on the cutout and raw data don't look matched.

What I would suggest is that you could set a breakpoint in tapas_physio_read_physlogfiles_siemens at about line 89 (hasScanTimingDicomImage) and step through from there to see whether the dicomHeader and its volume start/stop timestamps are read correctly and whether the timing information (until line 131) from scan and physiological logfile all seem consistent.

If that does not help, you might send me the data for this batch as well (.dcm, .puls, .resp file) to try it out quickly myself.

All the best, Lars

Dear Lars,

Above I told that the PMU logging was recorded for 5 runs. That is not the case. I started recordings just before my two functional T2*-scans which were placed at the end of my scanning session.

Thank you! Alex

Dear Alex,

What I would suggest is that you could set a breakpoint in tapas_physio_read_physlogfiles_siemens at about line 89 (hasScanTimingDicomImage) and step through from there to see whether the dicomHeader and its volume start/stop timestamps are read correctly and whether the timing information (until line 131) from scan and physiological logfile all seem consistent.

If that does not help, you might send me the data for this batch as well (.dcm, .puls, .resp file) to try it out quickly myself.

All the best, Lars

Dear Lars,

I ran tapas_physio_read_physlogfiles_siemens section-wise. Giving it the

`log_files = 

    respiration: 'E:\rm3\PMU\20191023_V02.resp'
        cardiac: 'E:\rm3\PMU\20191023_V02.puls'
    scan_timing: E:\rm3\V02\20191023_09h30_0011034185_dicom\0038061729\1.3.12.2.1107.5.2.43.67036.2019102310135153063141733.0.0.0\MR.27842.1611'`

the script runs correctly from line 89 – 104 and reads the dicom header via the function dicomHeader

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

It returns the appropriate value of 37412,4125000000 as read from the dicom tag for AquisitionTime.

Giving the script

log_files.align_scan ='last'
explicit_relative_start_acquisition = 0
cardiac_modality = 'PPU'
dt = []

the next vital section (lines 111 – 159) for cardiac read in runs and returns

logFooter. LogStartTimeSeconds 3.514548500000000e+04
logFooter.LogStopTimeSeconds   3.742517700000000e+04
logFooter.ScanStartTimeSeconds 35144832
logFooter.ScanStopTimeSeconds  37425315

which – after coversion to seconds - seems reasonable as the physiological recordings startet before and ended after MRI recordings. Thus, important requirements for syncronization are available. Nevertheless, the cutout in the physiological recordings fort eh fMRI runs is not in the correct position. Other good things appear in workspace correctly as far as I can see:

Now anyways, the whole physio script, as started from the batch editor, runs completely (as it did from the beginning) and returns

Running 'TAPAS PhysIO Toolbox'
**No scan trigger events provided**
    maxscan (incl. dummies) = 384 
    tmin (1st scan start (1st dummy))= 1721.14 s
    tmin (1st scan start (after dummies))= 1724.14 s
    tmax = 2297.07 s 
    mean TR =   1.50 s

Done    'TAPAS PhysIO Toolbox'
Done

in the command window.

The statement „No scan trigger events provided“ is derived from the script tapasphysioplot_raw_physdata.m from the folder …\PhysIO\code\sync. Nevertheless, the syncronization is not just a wrong plot as the resulting GLM regressors are the same for different onset-scan-triggers. At this point I guess that ons_secs.t is somehow not correct because the lines 51 – 56 in tapas_physio_plot_raw_physdata.m result in the else-statement:

`if has_scan_triggers
        plot(ons_secs.t, amp*ons_secs.acq_codes/max(ons_secs.acq_codes), 'c'); hold all;
        lg{end+1} = 'Scan trigger events';
    else
        verbose = tapas_physio_log('No scan trigger events provided', verbose, 0);
end
`

So, are there any other things to try to get the syncronization to be correct?

Thank you, Alex

Dear Lars,

to rule out that the incorrect syncronization is only due to the SIEMENS .resp and .puls - recordings, I got another data set with BIOPAC-recordings. Both data sets were collected at the same scanner. I was surprised to see the same behavior of the toolbox s described above in posting #3. The cutout region is at the very beginning or at the very end of the physiological recordings depending on matlabbatch{1}.spm.tools.physio.log_files.align_scan is set first or last. Of course, the scans of interest are in the middle of the physiological recordings. Here the batch script?

%-----------------------------------------------------------------------
% Job saved on 18-Dec-2019 14:59:11 by cfg_util (rev $Rev: 6460 $)
% spm SPM - SPM12 (6470)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.tools.physio.save_dir = {''};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Biopac_Mat';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'D:\Studie193_PLAU\Auswertung\Athletes Brain_Physiodaten_MRT\KP\KP01_MAAN88.mat'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'D:\Studie193_PLAU\Auswertung\Athletes Brain_Physiodaten_MRT\KP\KP01_MAAN88.mat'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'D:\Studie193_PLAU\Rohdaten\KP\KP01\04_Run1_Calibrierung\MR.19039.2487'};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 30;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 3;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 205;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 10;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'ECG';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

Considering this, the syncronization problem must have a more general cause, independent of the vendor type matlabbatch{1}.spm.tools.physio.log_files.vendor , in my case.

Do you know any help?

All the best! Alex

Dear Lars,

after trying different machines, different Tapas releases and different datasets with different physioloigical logging methods, I have one more thought to illustrate the still remaining one and the same problem:

In my example I have one physiological recording covering two fMRI runs. The first run covers the volumes MR.27842.264 to MR.27842.1229, and the second run MR.27842.1230 to MR.27842.1611.
To keep things simple and compareable, I set matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 100; in both calculations.

For run1 the batch looks like that:

matlabbatch{1}.spm.tools.physio.save_dir = {'D:\rm3\V02\physio_whole_head'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Siemens';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'D:\rm3\PMU\20191023_V02.puls'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'D:\rm3\PMU\20191023_V02.resp'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'D:\rm3\V02\20191023_09h30_0011034185_dicom\0038061729\1.3.12.2.1107.5.2.43.67036.2019102309451344116504598.0.0.0\MR.27842.1229'};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 12;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 1.5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 100;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = 0.125;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.scan_timing_log = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.max_heart_rate_bpm = 90;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

for run2 the only change is made in the dicom scan, i. e. the last scan volume of run2, which is MR.27842.1611 matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'D:\rm3\V02\20191023_09h30_0011034185_dicom\0038061729\1.3.12.2.1107.5.2.43.67036.2019102309451344116504598.0.0.0\MR.27842.1611'};

For run1 the command window says

Running job #1

Running 'TAPAS PhysIO Toolbox' No scan trigger events provided maxscan (incl. dummies) = 100 tmin (1st scan start (1st dummy))= 2847.15 s tmin (1st scan start (after dummies))= 2847.15 s tmax = 2997.03 s mean TR = 1.50 s

Note: Guessed 1 additional cardiac pulse(s) at time series end for phase estimation Done 'TAPAS PhysIO Toolbox' Done

For run2 the command window tells

Running job #1

Running 'TAPAS PhysIO Toolbox' No scan trigger events provided maxscan (incl. dummies) = 100 tmin (1st scan start (1st dummy))= 2147.14 s tmin (1st scan start (after dummies))= 2147.14 s tmax = 2297.01 s mean TR = 1.50 s

Done 'TAPAS PhysIO Toolbox' Done

Incorrectly,

tmin (1st scan start (1st dummy))= 2847.15 s

of run1 is later in timen than

tmin (1st scan start (1st dummy))= 2147.14 s

of run 2.

That results in these cutout regions: run1: run1_co

run2: run2_co

with RETROICOR regressors for run1: run1_RETROICOR

and run2: run2_RETROICOR

which are similar, but not the same due to a negligible estimation error, I guess.

As far as I can see, the cutout region does always include the very last data points of the recording, when matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last'; although physiological recordings were stopped some time after the last fmri run. Moreover, the first run ended substantially earlier than run2. Nevertheless, the cutout region is drawn at the very end of the physiological recordings.

When matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first'; the cutout region does always include the very first data points of the recording although fMRI run1 started some time after physiological recordings were switched on.

Hope that makes the syncronization failure a bit more clear. As mentioned above, this problem is the same with different data sets and tools. Just the scanner SIEMENS magnetom Prisma was the same for all data collection.

All the best, Alex

Dear Alex,

thank you so much for your patience and perseverance in solving this issue with me. I have now checked our bugfix, added some integration tests for your example data and ran all existing tests to make sure we broke nothing else :-)

The bottom line is: Indeed, there was a bug including DICOM timing information for Siemens VB logfiles, which was exacerbated for long physiological logfiles that lasted several runs or the whole scan session. Thanks to your example data Siemens_VB/PPU3T_Sync_First and Siemens_VB/PPU3T_Sync_Last, this was now resolved and is available as v7.2.2 release of PhysIO.

For now, this version can be found in the development branch of TAPAS, and will be included into the next release. Therein, also the new example data will be published.

All the best, Lars

translationalneuromodeling / tapas

SIEMENS PMU logging synchronization #55

Running job #1

Running job #1