Reviewer #3: (Daniele Marinazzo)

Dear authors,

Thanks for the opportunity to read this nice paper. I appreciated the quality of the data analysis, and the quest towards associating electrophysiology and BOLD data through a data-driven transfer function, which can be interpreted as a proxy of the HRF. Also I completely agree with you that we need to move beyond a canonical response.

There are a few issues I would like to discuss with you. I have done quite some work in this sense. On one hand this is good (and I think it's also the reason why I was invited to review this paper), on the other one there is always the risk that I have shaped my own goggles in these last years, and that I am projecting on your work some doubts and issues that I have with my own. In this case I apologize in advance, and I hope that we can have an enriching conversation.

Please forgive me if I start by my own work; there is always the danger that reviewers try to make authors write the paper that they would write themselves, I will keep this in mind, but on the other hand I think that the best way to convey my thoughts to you is to let them flow as they come.

So, here's our toolbox: https://www.nitrc.org/projects/rshrf. The idea behind it is that we can take peaks in the BOLD signal and take them as signatures of a pseudo neural event happening some time before at the neural level. This is in line with this work (which could also be relevant with respect to your power law figures):

Tagliazucchi E, Balenzuela P, Fraiman D, Chialvo DR. Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. Front Physiol. 2012;3:15. Published 2012 Feb 8. doi:10.3389/fphys.2012.00015 and with the subsequent spatial clustering approach which has been called coactivation patterns (CAP)

Liu X, Zhang N, Chang C, Duyn JH. Co-activation patterns in resting-state fMRI signals. Neuroimage. 2018;180(Pt B):485-494. doi:10.1016/j.neuroimage.2018.01.041 and innovation CAPs

Karahanoğlu FI, Caballero-Gaudes C, Lazeyras F, Van de Ville D. Total activation: fMRI deconvolution through spatio-temporal regularization. Neuroimage. 2013;73:121-134. doi:10.1016/j.neuroimage.2013.01.067 Karahanoğlu FI, Van De Ville D. Transient brain activity disentangles fMRI resting-state dynamics in terms of spatially and temporally overlapping networks. Nat Commun. 2015;6:7751. Published 2015 Jul 16. doi:10.1038/ncomms8751

Zoller DM, Bolton TAW, Karahanoglu FI, Eliez S, Schaer M, Van De Ville D. Robust Recovery of Temporal Overlap Between Network Activity Using Transient-Informed Spatio-Temporal Regression. IEEE Trans Med Imaging. 2019;38(1):291-302. doi:10.1109/TMI.2018.2863944

We then fit these peaks with a GLM, with the time lag as a free parameter. We use several families of basis functions. In the original paper (Wu GR, Liao W, Stramaglia S, Ding JR, Chen H, Marinazzo D. A blind deconvolution approach to recover effective connectivity brain networks from resting state fMRI data. Med Image Anal. 2013;17(3):365-374. doi:10.1016/j.media.2013.01.003) we used canonical HRF and FIR (together with the rBETA, which is basically the portion of the BOLD peak exceeding a certain threshold, as in the Tagliazucchi paper above).

We then included a mixture of gamma functions together with other families of basis functions in subsequent versions of the toolbox. Then we set up for validation of the approach with electrophysiological signatures, and that's where the doubts and pain kicked in. Some results on simultaneous EEG-fMRI, reported here (Wu G, Marinazzo D. 2015. Retrieving the Hemodynamic Response Function in resting state fMRI: methodology and applications. PeerJ PrePrints 3:e1317v1 https://doi.org/10.7287/peerj.preprints.1317v1 Wu GR, Deshpande G, Laureys S, Marinazzo D. Retrieving the Hemodynamic Response Function in resting state fMRI: Methodology and application. Conf Proc IEEE Eng Med Biol Soc. 2015;2015:6050-6053. doi:10.1109/EMBC.2015.7319771) were encouraging: for example we saw that the positive correlation between envelope of EEG and BOLD in the occipital cortex becomes more positive when we use instead the deconvolved BOLD and the EEG, while the negative correlation in the thalamus becomes more negative.

Other things present in the PeerJ preprint were encouraging too (and I mention them since I think that they can be relevant to the validation of your approach): namely the retrieval of a simulated ground truth HRF within certain realistic limits of SNR and jitter, the correlation with cerebral blood flow (even though physiological regressors should always be taken into account, see: Wu GR, Marinazzo D. Sensitivity of the resting-state haemodynamic response function estimation to autonomic nervous system fluctuations. Philos Trans A Math Phys Eng Sci. 2016;374(2067):20150190. doi:10.1098/rsta.2015.0190 and this becomes even more relevant when considering aging and clinical datasets), and some similarity across resting state networks.

So, the question is: can we really trust that peaks in M/EEG reflect the local pseudo-events that would originate the BOLD signal? Reading work by people who had thoroughly investigated this, e.g.

Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A. Neurophysiological investigation of the basis of the fMRI signal. Nature. 2001;412(6843):150-157. doi:10.1038/35084005

Chen X, Sobczak F, Chen Y, et al. Mapping optogenetically-driven single-vessel fMRI with concurrent neuronal calcium recordings in the rat hippocampus. Nat Commun. 2019;10(1):5239. Published 2019 Nov 20. doi:10.1038/s41467-019-12850-x

Yu X, He Y, Wang M, et al. Sensory and optogenetically driven single-vessel fMRI. Nat Methods. 2016;13(4):337-340. doi:10.1038/nmeth.3765

and conversing with them, I got (almost) convinced that it's unlikely that spikes in coarsely recorded or reconstructed M/EEG signal can be one to one mapped to the HRF inducing events that we use in GLM (calcium or even better glutamate signal could be a better choice).

Now, I like the way you associated HMM states with hemodynamic ones, thus adopting a more systemic/dynamical view, and taking fractional occupancy as a trigger. Do you think that these triggers can be better markers of BOLD-inducing neural events?

Other issues:

• What to make of events that are very close, and that would thus violate the assumption of linearity of the GLM?

• Apart from hemodynamic changes, can aging be associated with different electrophysiological spectral features (both periodic and aperiodic), which in turn could influence the HMM analysis?

• Detection of brain-behavior relationships with a non-huge dataset can be misleading, see for example this recent study:

Towards Reproducible Brain-Wide Association Studies Scott Marek, Brenden Tervo-Clemmens, Finnegan J. Calabro, David F. Montez, Benjamin P. Kay, Alexander S. Hatoum, Meghan Rose Donohue, William Foran, Ryland L. Miller, Eric Feczko, Oscar Miranda-Dominguez, Alice M. Graham, Eric A. Earl, Anders J. Perrone, Michaela Cordova, Olivia Doyle, Lucille A. Moore, Greg Conan, Johnny Uriarte, Kathy Snider, Angela Tam, Jianzhong Chen, Dillan J. Newbold, Annie Zheng, Nicole A. Seider, Andrew N. Van, Timothy O. Laumann, Wesley K. Thompson, Deanna J. Greene, Steven E. Petersen, Thomas E. Nichols, B.T. Thomas Yeo, Deanna M. Barch, Hugh Garavan, Beatriz Luna, Damien A. Fair, Nico U.F. Dosenbach bioRxiv 2020.08.21.257758; doi: 10.1101/2020.08.21.257758

• Why the parcellation in 38 regions? How are the results consistent/robust with finer parcellations?

• You state that the DMN "is susceptible to aging and neurodegenerative disease". That's certainly probable, the thing is that DMN is possibly sensitive to everything and specific to a very few things.

• Instead of a point-by-point statistical test, you could use the 3dMVM algorithm in AFNI (your reference 20) to test differences in the shape as a whole.

• You analyse data from older subjects only. How confident can you be that you are observing effects specific to aging?

Thanks for listening to this review version of "more of a comment than a question".