neuromodes.waves.simulate_waves

neuromodes.waves.simulate_waves(emodes: ArrayLike, evals: ArrayLike, nt: int = 1000, bold_out: bool = False, ext_input: ArrayLike | None = None, dt: float = 0.0001, r: float = 18.0, gamma: float = 116.0, pde_method: str = 'fourier', decomp_method: str = 'project', mass: spmatrix | ArrayLike | None = None, speed_limits: tuple[float, float] | None = (0, 150), scaled_hetero: ArrayLike | None = None, check_ortho: bool = True, seed: int | None = None, cache_input: bool = False, **balloon_params) → NDArray

Simulate neural activity or BOLD signals on the surface mesh using the NFT wave equation [1][2] using homogeneous or heterogeneous modes [3]. Optionally, the output of the wave model can be transformed into BOLD signals using the Balloon-Windkessel model [4][5].

Parameters

emodesarray-like

The eigenmodes array of shape (n_verts, n_modes), where n_verts is the number of vertices and n_modes is the number of eigenmodes.

evalsarray-like

The eigenvalues array of shape (n_modes,).

ntint, optional

Number of time points to simulate Default is 1000.

bold_outbool, optional

If True, simulate BOLD signal using the balloon model. If False, simulate neural activity. Default is False.

ext_inputarray-like, optional

External input array of shape (n_verts, n_timepoints). If None, random input is generated. Default is None.

dtfloat, optional

Time step for simulation in seconds. Default is 1e-4.

rfloat, optional

Spatial length scale of wave propagation in millimeters. Default is 18.0.

gammafloat, optional

Damping rate of wave propagation in seconds^-1. Default is 116.0.

pde_methodstr, optional

Method for solving the wave PDE. Either ‘fourier’ or ‘ode’. Default is ‘fourier’.

decomp_methodstr, optional

The method used for the eigendecomposition, either ‘project’ to project data into a mass-orthonormal space or ‘regress’ for least-squares fitting. Note that the beta values from ‘regress’ tend towards those from ‘project’ when more modes are provided. Default is ‘project’.

massarray-like, optional

The mass matrix of shape (n_verts, n_verts) used for the decomposition when method is ‘project’. If using EigenSolver, provide its self.mass. Default is None.

speed_limitstuple, optional

If any wave speeds are outside this range (in m/s), a warning is raised. If None, no warning is raised. Default is (0, 150).

scaled_heteroarray-like, optional

Scaled heterogeneity map of shape (n_verts,), used only to check wave speeds (see speed_limits above). If not provided, wave speed is assumed to be spatially uniform. To scale a heterogeneity map, use the eigen.scale_hetero function. Default is None.

check_orthobool, optional

Whether to check if emodes are mass-orthonormal before using the ‘project’ method for decomposition. Default is True.

seedint, optional

Random seed for generating external input. Default is None.

cache_inputbool, optional

If True and ext_input is None, cache the generated random input to avoid recomputation for the same values of nt, seed, and number of rows (vertices) in emodes. Inputs are cached in the directory specified by the CACHE_DIR environment variable. If not set, the user’s home directory is chosen. Default is False.

**balloon_params

Optional balloon model parameters to override defaults (e.g., rho, k1). See get_balloon_params() for available parameters. Only used when bold_out=True.

Returns

np.ndarray

Simulated neural activity or BOLD signal of shape (n_verts, n_timepoints).

Raises

ValueError

If emodes does not have shape (n_verts, n_modes), where n_verts ≥ n_modes.

ValueError

If evals does not have shape (n_modes,).

ValueError

If r, gamma, or dt is not positive.

ValueError

If nt is not a positive integer.

ValueError

If speed_limits is not a tuple (min_speed, max_speed), where 0 ≤ min_speed < max_speed.

ValueError

If ext_input is provided but does not have shape (n_verts, nt).

ValueError

If pde_method is not ‘fourier’ or ‘ode’.

RuntimeError

If the ODE solver fails when using pde_method=’ode’ and bold_out=True.

Notes

Since the simulation begins at rest, consider discarding the first ~50 timepoints to allow the system to reach a steady state.

References

[1]

Robinson, P. A., et al. (1997). Propagation and stability of waves of electrical activity in the cerebral cortex. Physical Review E. https://doi.org/10.1103/physreve.56.826

[2]

Pang, J. C., et al. (2023). Geometric constraints on human brain function. Nature. https://doi.org/10.1038/s41586-023-06098-1

[3]

Barnes, V., et al. (2026). Regional heterogeneity shapes macroscopic wave dynamics of the human and non-human primate cortex. bioRxiv. https://doi.org/10.64898/2026.01.22.701178

[4]

Buxton, R. B., et al. (1998). Dynamics of blood flow and oxygenation changes during brain activation: The balloon model. Magnetic Resonance in Med. https://doi.org/10.1002/mrm.1910390602

[5]

Stephan, K. E., et al. (2007). Comparing hemodynamic models with DCM. NeuroImage. https://doi.org/10.1016/j.neuroimage.2007.07.040