OSCR

Estimating fMRI timescale maps.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Methods › Estimation of standard errors › Time-domain standard error estimator ↔ fmri_timescales/timescale_utils.py, lines 133–256 · score 0.64 · Bartlett kernel, lag truncation, weighted, regression, standard errors, domain
  2. [2] § Simulations › Simulation settings ↔ fmri_timescales/sim.py, lines 9–62 · score 0.55 · multivariate normal distribution, Cholesky, matrix, Toeplitz, simulate, ACFs

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 256 lines · 9.2 KB · MIT · 1 match

  1. from typing import Optional
  3. import numpy as np
  4. from joblib import Parallel, delayed
  5. from scipy.optimize import curve_fit
  6. from sklearn.base import BaseEstimator
  8. from fmri_timescales import acf_utils
  11. def newey_west_omega(u: np.ndarray, n_lags: int | None = None) -> float:
  12. n_u = len(u)
  13. if n_lags is None:
  14. n_lags = int(np.floor(4 * (n_u / 100.0) ** (2 / 9)))
  15. weights = 1 - np.arange(n_lags + 1) / (n_lags + 1)
  16. omega = weights[0] * np.sum(u**2)
  17. for lag in range(1, n_lags + 1):
  18. omega += weights[lag] * (2 * np.sum(u[lag:] * u[:-lag]))
  19. return omega
  22. def _phi_to_tau(phis: np.ndarray, se_phis: np.ndarray, lag: int = 1) -> tuple:
  23. """phi to tau (timescale), and apply delta method to std err"""
  24. phis_abs = np.abs(phis) # tau undefined for negative phi
  25. taus = -lag / np.log(phis_abs)
  26. se_taus = (lag / (phis_abs * np.log(phis_abs) ** 2)) * se_phis
  27. return taus, se_taus
  30. class TD(BaseEstimator):
  31. """Time Domain (TD) Linear Model, Fit by Linear Least Squares.
  33. Parameters
  34. ----------
  35. var_estimator : str, optional
  36. The variance estimator to use. Options are "newey-west" or "non-robust", by default "newey-west"
  37. var_n_lags : int, optional
  38. The lag truncation number for the bartlett kernel, by default None
  39. copy_X : bool, optional
  40. If True X will be copied, else it may be overwritten, by default False
  41. n_jobs : int, optional
  42. The number of jobs to use for the computation, by default None
  44. Attributes
  45. ----------
  46. estimates_ : dict
  47. A dictionary containing four np.ndarray of shape (n_regions, ):
  48. - "phi": AR(1) coefficient estimates, for each region in X.
  49. - "se(phi)": Standard errors of AR(1) coefficients.
  50. - "tau": Timescale estimates, for each region in X.
  51. - "se(tau)": Standard errors of timescales.
  53. Examples
  54. --------
  55. >>> from fmri_timescales import sim, timescale_utils
  56. >>> X = sim.sim_ar(ar_coeffs=[0.8], n_timepoints=1000) # x_t = 0.8 x_{t-1} + e_t
  57. >>> td = timescale_utils.TD(var_estimator="newey-west", var_n_lags=10)
  58. >>> td.fit(X=X, n_timepoints=1000).estimates_
  59. {'phi': array([0.79789847]), 'se(phi)': array([0.02045074]), 'tau': array([4.42920958]), 'se(tau)': array([0.50282146])}
  60. """
  62. def __init__(
  63. self,
  64. var_estimator: str = "newey-west",
  65. var_n_lags: int | None = None,
  66. copy_X: bool = False,
  67. n_jobs: int | None = None,
  68. ) -> None:
  69. self.var_estimator = var_estimator
  70. self.var_n_lags = var_n_lags
  71. self.copy_X = copy_X
  72. self.n_jobs = n_jobs
  74. @delayed
  75. def _fit_td(self, x: np.ndarray) -> tuple:
  76. """fit model to a single timeseries x in X"""
  77. T = len(x) - 1
  78. # x_t = X[1:], x_{t-1} = x[:-1] (Hz=1)
  79. phi_ = np.sum(x[1:] * x[:-1]) / np.sum(x[:-1] ** 2)
  81. # variance estimators
  82. def non_robust():
  83. e_ = x[1:] - phi_ * x[:-1]
  84. q_ = np.sum(x[:-1] ** 2)
  85. sigma2_ = (1 / T) * np.sum(e_**2)
  86. return (1 / q_) * sigma2_
  88. def newey_west():
  89. e_ = x[1:] - phi_ * x[:-1]
  90. q_ = np.sum(x[:-1] ** 2)
  91. u_ = x[:-1] * e_
  92. omega_ = newey_west_omega(u_, n_lags=self.var_n_lags)
  93. return (1 / q_) * omega_ * (1 / q_)
  95. var_estimators = {"non-robust": non_robust, "newey-west": newey_west}
  97. if self.var_estimator not in var_estimators:
  98. raise ValueError("var_estimator must be either 'newey-west' or 'non-robust'")
  100. var_ = var_estimators[self.var_estimator]()
  101. return phi_, np.sqrt(var_)
  103. def fit(self, X: np.ndarray, n_timepoints: int):
  104. """Fit the TD model.
  106. Parameters
  107. ----------
  108. X : np.ndarray of shape (n_timepoints, n_regions)
  109. An array containing the timeseries of each region.
  110. n_timepoints : int
  111. The number of timepoints in X.
  113. Raises
  114. ------
  115. ValueError
  116. If `X` is not in (n_timepoints, n_regions) form.
  117. """
  119. if X.ndim != 2 or X.shape[0] != n_timepoints:
  120. raise ValueError("X should be in (n_timepoints, n_regions) form")
  121. X = X.copy() if self.copy_X else X
  122. X = (X - X.mean(axis=0)) / X.std(axis=0) # mean zero, variance 1
  124. with Parallel(n_jobs=self.n_jobs) as parallel:
  125. td_fits = parallel(self._fit_td(X[:, idx]) for idx in range(X.shape[1]))
  126. phis_, se_phis_ = map(np.array, zip(*td_fits))
  127. taus_, se_taus_ = _phi_to_tau(phis_, se_phis_)
  129. self.estimates_ = {"phi": phis_, "se(phi)": se_phis_, "tau": taus_, "se(tau)": se_taus_}
  130. return self
  133. class AD(BaseEstimator):
  134. """Autocorrelation Domain (AD) Nonlinear Model, fit by Nonlinear Least Squares.
  136. Parameters
  137. ----------
  138. var_estimator : str, optional
  139. The variance estimator to use. Options are "newey-west" or "non-robust", by default "newey-west"
  140. var_n_lags : int, optional
  141. The lag truncation number for the bartlett kernel, by default None
  142. acf_n_lags : int, optional
  143. The lag truncation number for the autocorrelation function, by default None
  144. copy_X : bool, optional
  145. If True X will be copied, else it may be overwritten, by default False
  146. n_jobs : _type_, optional
  147. The number of jobs to use for the computation, by default None
  149. Attributes
  150. ----------
  151. estimates_ : dict
  152. A dictionary containing two np.ndarray of shape (n_regions, ):
  153. - "tau": Timescale estimates, for each region in X.
  154. - "se(tau)": Standard errors of timescales.
  156. Examples
  157. --------
  158. >>> from fmri_timescales import sim, timescale_utils
  159. >>> X = sim.sim_ar(ar_coeffs=[0.8], n_timepoints=1000) # x_t = 0.8 x_{t-1} + e_t
  160. >>> ad = timescale_utils.AD(var_estimator="newey-west", var_n_lags=10, acf_n_lags=50)
  161. >>> ad.fit(X=X, n_timepoints=1000).estimates_
  162. {'phi': array([0.78021651]), 'se(phi)': array([0.02814532]), 'tau': array([4.02927146]), 'se(tau)': array([0.58565806])}
  163. """
  165. def __init__(
  166. self,
  167. var_estimator: str = "newey-west",
  168. var_n_lags: int | None = None,
  169. acf_n_lags: int | None = None,
  170. copy_X: bool = False,
  171. n_jobs: int | None = None,
  172. ) -> None:
  173. self.var_estimator = var_estimator
  174. self.var_n_lags = var_n_lags
  175. self.acf_n_lags = acf_n_lags
  176. self.copy_X = copy_X
  177. self.n_jobs = n_jobs
  179. @delayed
  180. def _fit_ad(self, x: np.ndarray) -> tuple:
  181. """fit model to a single timeseries/autocorrelation function x in X"""
  183. T = len(x)
  185. # acf estimator
  186. x_acf = acf_utils.ACF(n_lags=self.acf_n_lags + 1).fit_transform(x.reshape(-1, 1), T).squeeze()[1:]
  188. # regression function (m), and its linearized regressor (dm_dphi)
  189. ks = np.arange(1, len(x_acf) + 1)
  191. def m(ks, phi):
  192. return phi**ks
  194. def jac(ks, phi):
  195. return (ks * phi ** (ks - 1)).reshape(-1, 1)
  197. # phi estimator
  198. eps = 1e-10
  199. phi_, _ = curve_fit(f=m, xdata=ks, ydata=x_acf, p0=1e-2, bounds=(-1 + eps, +1 - eps), ftol=1e-6, jac=jac)
  200. phi_ = phi_.squeeze()
  202. # variance estimators
  203. def non_robust():
  204. e_ = x[1:] - phi_ * x[:-1]
  205. q_ = np.sum(x[:-1] ** 2)
  206. sigma2_ = (1 / T) * np.sum(e_**2)
  207. return (1 / q_) * sigma2_
  209. def newey_west():
  210. q_ = np.sum((ks * phi_ ** (ks - 1)) ** 2)
  211. weights = ks * (phi_ ** (ks - 1))
  212. weights_phi = weights * (phi_**ks)
  213. conv1 = np.convolve(x, weights, mode="full")
  214. conv2 = np.convolve(x**2, weights_phi, mode="full")
  215. u_ = x[len(ks) : T] * conv1[len(ks) - 1 : T - 1] - conv2[len(ks) - 1 : T - 1]
  216. omega_ = (1 / len(u_) ** 2) * newey_west_omega(u_, n_lags=self.var_n_lags)
  217. return (1 / q_) * omega_ * (1 / q_)
  219. var_estimators = {"non-robust": non_robust, "newey-west": newey_west}
  221. if self.var_estimator not in var_estimators:
  222. raise ValueError("var_estimator must be either 'newey-west' or 'non-robust'")
  224. var_ = var_estimators[self.var_estimator]()
  225. return phi_, np.sqrt(var_)
  227. def fit(self, X: np.ndarray, n_timepoints: int):
  228. """Fit the AD model.
  230. Parameters
  231. ----------
  232. X : np.ndarray of shape (n_timepoints, n_regions)
  233. An array containing the timeseries of each region.
  234. n_timepoints : int
  235. The number of timepoints in X.
  237. Raises
  238. ------
  239. ValueError
  240. If `X` is not in (n_timepoints, n_regions) form.
  241. """
  243. if X.ndim != 2 or X.shape[0] != n_timepoints:
  244. raise ValueError("X should be in (n_timepoints, n_regions) form")
  245. X = X.copy() if self.copy_X else X
  246. X = (X - X.mean(axis=0)) / X.std(axis=0) # mean zero, variance 1
  247. if self.acf_n_lags is None:
  248. self.acf_n_lags = n_timepoints // 100
  250. with Parallel(n_jobs=self.n_jobs) as parallel:
  251. ad_fits = parallel(self._fit_ad(X[:, idx]) for idx in range(X.shape[1]))
  252. phis_, se_phis_ = map(np.array, zip(*ad_fits))
  253. taus_, se_taus_ = _phi_to_tau(phis_, se_phis_)
  255. self.estimates_ = {"phi": phis_, "se(phi)": se_phis_, "tau": taus_, "se(tau)": se_taus_}
  256. return self

timescale_utils.py at commit cd0072e, under MIT · at the source

Overview

Authors: Gabriel Riegner1, Samuel Davenport2, Bradley Voytek1,3,4, Armin Schwartzman1,2
  1. Halicioğlu Data Science Institute, University of California San Diego, La Jolla, CA, United States
  2. Division of Biostatistics, University of California San Diego, La Jolla, CA, United States
  3. Department of Cognitive Science, University of California San Diego, La Jolla, CA, United States
  4. Neurosciences Graduate Program, University of California San Diego, La Jolla, CA, United States
Institutions: University of California San Diego (United States)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1248
Dates: received 23 April 2025; accepted 27 April 2026; published online 4 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1248 · PMID 42253608 · PMCID PMC13237991 · OpenAlex W4409838608
Open access: diamond, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.04.23.650300
Status: code verified
Categories: fMRI (modality), human (organism)
Methods: Connectivity, Smoothing, state filtering, decompositions, Statistics, fMRI & imaging, Spectral & time-frequency
Keywords: time-domain linear model, autocorrelation-domain nonlinear model, uncertainty quantification, statistical inference, human connectome project, functional brain organization
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: NIH (R01MH128923); Halicioglu Data Science Institute Fellowship
Citations: not cited yet (Europe PMC); 56 references in the paper

Abstract

Brain activity unfolds over hierarchical timescales that reflect how brain regions integrate and process information, linking functional and structural organization. While timescale studies are prevalent, existing estimation methods rely on the restrictive assumption of exponentially decaying temporal autocorrelation and only provide point estimates without standard errors, limiting statistical inference. In this paper, we formalize and evaluate two methods for mapping timescales in resting-state fMRI: a time-domain fit of an autoregressive (AR1) model and an autocorrelation-domain fit of an exponential decay model. Rather than assuming exponential autocorrelation decay, we define timescales by projecting the fMRI time series onto these approximating models, requiring only stationarity and mixing conditions while incorporating robust standard errors to account for model misspecification. We introduce theoretical properties of timescale estimators and show parameter recovery in realistic simulations, as well as applications to fMRI from the Human Connectome Project. Comparatively, the time-domain method produces more accurate estimates under model misspecification, remains computationally efficient for high-dimensional fMRI data, and yields maps aligned with known functional brain organization. In this work, we show valid statistical inference on fMRI timescale maps, and provide Python implementations of all methods.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

griegner/fmri-timescales

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: cd0072ee3b42f271e49de3682d5f7a41bb9070d2, 18 May 2026
Languages: Python (13), Jupyter (6)
Size: 115 files, 19 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, license file, environment (pyproject.toml), tests, continuous integration, documentation, 6 notebooks
Not found: CITATION.cff
Tools: NumPy (17 files), Matplotlib (8 files), SciPy (7 files), NiBabel (5 files), scikit-learn (4 files), statsmodels (4 files), Nilearn (3 files), neuromaps (1 file), pandas (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
21 files
At the source: github.com/griegner/fmri-timescales

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 19 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and Code Availability

All simulation results and fMRI timescale maps, inclusive of the code by which they were derived, can be accessed on https://github.com/griegner/fmri-timescales. The code is under the open source MIT license, allowing access and reuse with attribution. The Human Connectome Project young adult dataset (ages 22-35; 2018 release) used in this study is publicly accessible under a data usage agreement, which describes specific terms for data use and sharing.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 4 authors, 6 keywords, 2 funders, 55 references.

Cite

This paper

Riegner, G., Davenport, S., Voytek, B., & Schwartzman, A. (2026). Estimating fMRI timescale maps. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1248. https://doi.org/10.1162/imag.a.1248

BibTeX

@article{riegner2026estimating,
author = {Riegner, Gabriel and Davenport, Samuel and Voytek, Bradley and Schwartzman, Armin},
title = {{Estimating fMRI timescale maps}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = jun,
volume = {4},
pages = {IMAG.a.1248},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1248},
url = {https://doi.org/10.1162/imag.a.1248},
pmid = {42253608},
pmcid = {PMC13237991}
}

RIS

TY - JOUR
AU - Riegner, Gabriel
AU - Davenport, Samuel
AU - Voytek, Bradley
AU - Schwartzman, Armin
TI - Estimating fMRI timescale maps
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/06/04
VL - 4
SP - IMAG.a.1248
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1248
UR - https://doi.org/10.1162/imag.a.1248
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1248",
"type": "article-journal",
"title": "Estimating fMRI timescale maps",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Riegner",
"given": "Gabriel"
},
{
"family": "Davenport",
"given": "Samuel"
},
{
"family": "Voytek",
"given": "Bradley"
},
{
"family": "Schwartzman",
"given": "Armin"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1248",
"DOI": "10.1162/imag.a.1248",
"PMID": "42253608",
"PMCID": "PMC13237991",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1248",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
4
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-74466-2 [code]
Neuromorphic hierarchical modular reservoirs.
Journal: Nature communications
In common: neuromaps, Nilearn, NiBabel, 6 other tools, 11 references
[2] doi:10.1162/imag.a.1275 [code]
Glucose metabolism echoes long-range temporal correlations in the human brain.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: NiBabel, scikit-learn, pandas, 3 other tools, fMRI, 9 references
[3] doi:10.1038/s41540-026-00727-x [code]
Association-sensory spatiotemporal hierarchy and functional gradient-regularised recurrent neural network with implications for schizophrenia.
Journal: NPJ systems biology and applications
In common: statsmodels, pandas, SciPy, 2 other tools, 9 references
[4] doi:10.1093/cercor/bhag121 [code]
Preserved intrinsic neural timescale organization with hierarchical variation in autism spectrum disorder.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: 10 references
[5] doi:10.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: neuromaps, Nilearn, NiBabel, 6 other tools, 4 references
[6] doi:10.7554/elife.110294 [code]
Arousal modulates functional connectivity through structured and hemispherically asymmetric community architecture during wakefulness.
Journal: eLife
In common: neuromaps, NiBabel, statsmodels, 5 other tools, fMRI, 4 references
[7] doi:10.21203/rs.3.rs-9326213/v1 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: Research Square (preprint)
In common: neuromaps, Nilearn, NiBabel, 5 other tools, fMRI, 4 references
[8] doi:10.64898/2026.03.09.710558 [code]
Multi-task fMRI outperforms resting-state fMRI for revealing task-invariant organization of the human brain
Journal: bioRxiv (preprint)
In common: neuromaps, Nilearn, NiBabel, 5 other tools, fMRI, 4 references
[9] doi:10.1038/s41467-026-72931-6 [code]
Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.
Journal: Nature communications
In common: Nilearn, NiBabel, statsmodels, 5 other tools, fMRI, 4 references
[10] doi:10.1162/netn.a.547 [code]
An evaluation of the efficacy of single-echo and multi-echo fMRI denoising strategies.
Journal: Network neuroscience (Cambridge, Mass.)
In common: Nilearn, NiBabel, statsmodels, 5 other tools, fMRI, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.