API and documentation¶
C++ models¶
vbi.models.cpp.jansen_rit¶
- class JR_sde(par={})[source]¶
Jansen-Rit model C++ implementation.
- Parameters:
- par: dict
Including the following: - A : [mV] determine the maximum amplitude of the excitatory PSP (EPSP) - B : [mV] determine the maximum amplitude of the inhibitory PSP (IPSP) - a : [Hz] 1/tau_e, \(\sum\) of the reciprocal of the time constant of passive membrane and all other spatially distributed delays in the dendritic network - b : [Hz] 1/tau_i - r [mV] the steepness of the sigmoidal transformation. - v0 parameter of nonlinear sigmoid function - vmax parameter of nonlinear sigmoid function - C_i [list or np.array] average number of synaptic contacts in th inhibitory and excitatory feedback loops - noise_amp - noise_std
dt [second] integration time step
t_initial [s] initial time
t_end [s] final time
method [str] method of integration
t_transition [s] time to reach steady state
dim [int] dimention of the system
- valid_params = ['noise_seed', 'seed', 'G', 'weights', 'A', 'B', 'a', 'b', 'noise_mu', 'noise_std', 'vmax', 'v0', 'r', 'C0', 'C1', 'C2', 'C3', 'dt', 'method', 't_transition', 't_end', 'control', 'output', 'RECORD_AVG', 'initial_state']¶
- class JR_sdde(par={})[source]¶
- valid_params = ['weights', 'delays', 'dt', 't_end', 'G', 'A', 'a', 'B', 'b', 'mu', 'nstart', 't_end', 't_transition', 'sigma', 'C', 'record_step', 'C0', 'C1', 'C2', 'C3', 'vmax', 'r', 'v0', 'output', 'sti_ti', 'sti_duration', 'sti_amplitude', 'sti_gain', 'noise_seed', 'seed', 'method']¶
vbi.models.cpp.km¶
- class KM_sde(par)[source]¶
Kuramoto model with noise (sde), C++ implementation.
- Parameters:
- pardict
Dictionary of parameters.
- valid_parameters = ['G', 'dt', 'noise_amp', 'omega', 'weights', 'noise_seed', 'seed', 'alpha', 't_initial', 't_transition', 't_end', 'output', 'num_threads', 'initial_state', 'type']¶
vbi.models.cpp.mpr¶
- class MPR_sde(par: dict = {}, parbold={})[source]¶
MPR model
- run(par: dict = {}, x0: ndarray | None = None, verbose: bool = False)[source]¶
Integrate the MPR model with the given parameters.
- Parameters:
- pardict
Dictionary of parameters.
- x0array_like
Initial condition of the system.
- verbosebool
If True, print the progress of the simulation.
- Returns:
- boldarray_like
Simulated BOLD signal.
vbi.models.cpp.vep¶
vbi.models.cpp.wc¶
- class WC_ode(par={})[source]¶
References:
[WC_1972] (1,2)Wilson, H.R. and Cowan, J.D. Excitatory and inhibitory interactions in localized populations of model neurons, Biophysical journal, 12: 1-24, 1972.
[WC_1973]Wilson, H.R. and Cowan, J.D A Mathematical Theory of the Functional Dynamics of Cortical and Thalamic Nervous Tissue
[D_2011]Daffertshofer, A. and van Wijk, B. On the influence of amplitude on the connectivity between phases Frontiers in Neuroinformatics, July, 2011
Used Eqns 11 and 12 from [WC_1972] in
rhs. P and Q represent external inputs, which when exploring the phase portrait of the local model are set to constant values. However in the case of a full network, P and Q are the entry point to our long range and local couplings, that is, the activity from all other nodes is the external input to the local population [WC_1973], [D_2011] .The default parameters are taken from figure 4 of [WC_1972], pag. 10.
vbi.models.cpp.damp_oscillator¶
- class DO(par={})[source]¶
Damp Oscillator model class.
- valid_params = ['a', 'b', 'dt', 't_start', 't_end', 't_cut', 'initial_state', 'method', 'output']¶
- __init__(par={})[source]¶
- Parameters:
- pardictionary
parameters which includes the following: - dt [double] time step. - t_start [double] initial time for simulation. - t_end [double] final time for simulation. - initial_state [list] initial state of the system.
- update_par(par={})[source]¶
Update model parameters.
- Parameters:
- pardict, optional
Dictionary of parameters to update. Keys must be valid parameter names.
- run(par={}, x0=None, verbose=False)[source]¶
Run the damped oscillator simulation.
- Parameters:
- pardict, optional
Dictionary of parameters to update for this simulation run. Any parameter from the class documentation can be updated.
- x0array-like, optional
Initial state vector [x0, y0] of length 2. If None, uses the initial state set during initialization.
- verbosebool, optional
If True, print verbose output during simulation. Default is False.
- Returns:
- dict
Dictionary containing simulation results:
‘t’ : np.ndarray of shape (n_steps,) - time points
‘x’ : np.ndarray of shape (n_steps, 2) - simulated trajectories where x[:, 0] is x(t) and x[:, 1] is y(t)
Cupy models¶
vbi.models.cupy.mpr¶
- class MPR_sde(par: dict = {}, Bpar: dict = {})[source]¶
Montbrio-Pazo-Roxin model Cupy and Numpy implementation.
- Parameters:
- G: float. np.ndarray
global coupling strength
- dt: float
time step
- dt_bold: float
time step for Balloon model
- J: float, np.ndarray
model parameter
- eta: float, np.ndarray
model parameter
- tau:
model parameter
- delta:
model parameter
- tr: float
repetition time of fMRI
- noise_amp: float, np.array
amplitude of noise
- same_noise_per_sim:
same noise for all simulations
- iapp: float, np.ndarray
external input
- t_start: float
initial time
- t_cut: float
transition time
- t_end: float
end time
- num_nodes: int
number of nodes
- weights: np.ndarray
weighted connection matrix
- rv_decimate: int
sampling step from r and v
- output: str
output directory
- RECORD_TS: bool
store r and v time series
- num_sim: int
number of simulations
- method: str
integration method
- engine: str
cpu or gpu
- seed: int
seed for random number generator
- dtype: str
float or f
- initial_state: np.ndarray
initial state
- same_initial_state: bool
same initial state for all simulations
- set_initial_state(nn, ns, engine, seed=None, same_initial_state=False, dtype=<class 'float'>)[source]¶
Set initial state
- Parameters:
- nnint
number of nodes
- nsint
number of simulations
- enginestr
cpu or gpu
- same_initial_conditionbool
same initial condition for all simulations
- seedint
random seed
- dtypestr
float: float64 f : float32
vbi.models.cupy.ghb¶
vbi.models.cupy.jansen_rit¶
- class JR_sde(par: dict = {})[source]¶
Jansen-Rit model cupy implementation.
Parameters¶ Name
Explanation
Default Value
AExcitatory post synaptic potential amplitude.
3.25
BInhibitory post synaptic potential amplitude.
22.0
1/aTime constant of the excitatory postsynaptic potential.
a: 0.1 (1/a: 10.0)
1/bTime constant of the inhibitory postsynaptic potential.
b: 0.05 (1/b: 20.0)
C0,C1Average numbers of synapses between excitatory populations. If array-like, it should be the shape of (
num_nodes,num_sim).C0: 1.0 * 135.0, C1: 0.8 * 135.0
C2,C3Average numbers of synapses between inhibitory populations. If array-like, it should be the shape of (
num_nodes,num_sim).C2: 0.25 * 135.0, C3: 0.25 * 135.0
vmaxMaximum firing rate
0.005
vPotential at half of maximum firing rate
6.0
rSlope of sigmoid function at
v_00.56
GScaling the strength of network connections. If array-like, it should of length
num_sim.1.0
muMean of the noise
0.24
noise_ampAmplitude of the noise
0.01
weightsWeight matrix of shape (
num_nodes,num_nodes)None
num_simNumber of simulations
1
dtTime step
0.01
t_endEnd time of simulation
1000.0
t_cutCut time
500.0
engine“cpu” or “gpu”
“cpu”
method“heun” or “euler” method for integration
“heun”
seedRandom seed
None
initial_stateInitial state of the system of shape (
num_nodes,num_sim)None
same_initial_stateIf True, all simulations have the same initial state
False
same_noise_per_simIf True, all simulations have the same noise
False
decimateDecimation factor for the output time series
1
dtypeData type to use for the simulation,
floatforfloat64orfforfloat32.“float”
- check_parameters(par)[source]¶
Check if the provided parameters are valid.
- Parameters:
- pardict
Dictionary of parameters to check.
- Raises:
- ValueError
If any parameter in
paris not valid.
- get_default_parameters() dict[source]¶
Default parameters for the Jansen-Rit model
- Parameters:
- nnint
number of nodes
- Returns:
- paramsdict
default parameters
- S_(x)[source]¶
Compute the sigmoid function.
This function calculates the sigmoid of the input
xusing the parametersvmax,r, andv0.
- Parameters:
- xfloat or array-like
The input value(s) for which to compute the sigmoid function.
- Returns:
- float or array-like
The computed sigmoid value(s).
- f_sys(x0, t)[source]¶
Compute the derivatives of the Jansen-Rit neural mass model.
- Parameters:
- x0array_like
Initial state vector of the system. It should have a shape of (6*nn, ns), where nn is the number of neurons and ns is the number of simulations.
- tfloat
Current time point (not used in the computation but required for compatibility with ODE solvers).
- Returns:
- dxarray_like
Derivatives of the state vector. It has the same shape as
x0.- The function computes the derivatives of the state vector based on the Jansen-Rit model equations.
- euler(x0, t)[source]¶
Perform one step of the Euler-Maruyama method for stochastic differential equations.
- Parameters:
- x0array_like
The initial state of the system.
- tfloat
The current time.
- Returns:
- array_like
The updated state of the system after one Euler step.
- set_initial_state(nn, ns, engine, seed=None, same_initial_state=False, dtype='float')[source]¶
set initial state for the Jansen-Rit model
- Parameters:
- nn: int
number of nodes
- ns: int
number of simulations
- engine: str
cpu or gpu
- seed: int
random seed
- same_initial_state: bool
if True, all simulations have the same initial state
- dtype: str
data type
- Returns:
- y0: array [nn, ns]
initial state
vbi.models.cupy.utils¶
- get_module(engine='gpu')[source]¶
Switches the computational engine between GPU and CPU.
- Parameters:
- enginestr, optional
The computational engine to use. Can be either “gpu” or “cpu”. Default is “gpu”.
- Returns:
- module
The appropriate array module based on the specified engine. If “gpu”, returns the CuPy module. If “cpu”, returns the NumPy module.
- Raises:
- ValueError
If the specified engine is not “gpu” or “cpu”.
If CuPy is not installed.
- tohost(x)[source]¶
move data to cpu if it is on gpu
- Parameters:
- x: array
data
- Returns:
- array
data moved to cpu
- is_on_cpu(x)[source]¶
Check if input is on CPU (i.e., not a CuPy array)
- Parameters:
- x: any
Input to check
- Returns:
- bool
True if input is not a CuPy array (i.e., is on CPU), False otherwise
- is_on_gpu(x)[source]¶
Check if input is on GPU (i.e., is a CuPy array)
- Parameters:
- x: any
Input to check
- Returns:
- bool or None
True if input is a CuPy array (i.e., is on GPU) False if input is not a CuPy array None if CuPy is not installed
- repmat_vec(vec, ns, engine)[source]¶
repeat vector ns times
- Parameters:
- vec: array 1d
vector to be repeated
- ns: int
number of repetitions
- engine: str
cpu or gpu
- Returns:
- vec: array [len(vec), n_sim]
repeated vector
- is_seq(x)[source]¶
check if x is a sequence
- Parameters:
- x: any
variable to be checked
- Returns:
- bool
True if x is a sequence
- prepare_vec(x, ns, engine, dtype='float')[source]¶
check and prepare vector dimension and type
- Parameters:
- x: array 1d
vector to be prepared, if x is a scalar, only the type is changed
- ns: int
number of simulations
- engine: str
cpu or gpu
- dtype: str or numpy.dtype
data type to convert to
- Returns:
- x: array [len(x), n_sim]
prepared vector
- get_(x, engine='cpu', dtype='f')[source]¶
- Parameters:
- xarray-like
The input array to be converted.
- enginestr, optional
The computation engine to use. If “gpu”, the array is transferred from GPU to CPU. Defaults to “cpu”.
- dtypestr, optional
The desired data type for the output array. Defaults to “f”.
- Returns:
- array-like
The converted array with the specified data type.
- dtype_convert(dtype)[source]¶
Convert a string representation of a data type to a numpy dtype.
- Parameters:
- dtypestr
The string representation of the data type (e.g., “float”, “float32”, “float64”).
- Returns:
- numpy.dtype
The corresponding numpy dtype.
- Raises:
- ValueError
If the input string does not match any known data type.
vbi.models.cupy.bold¶
vbi.models.cupy.km¶
- class KM_sde(par={})[source]¶
Kuramoto model with noise (stochastic differential equation)
- Parameters:
- G: float
global coupling strength
- dt: float
time step
- noise_amp: float
noise amplitude
- weights: array
weighted connection matrix
- omega: array
natural angular frequency
- seed: int
fix random seed for initial state
- t_cut: float
transition time
- t_end: float
end time
- decimate: int
decimate the output time series
- output: str
output directory
- initial_state: array
initial state
- engine: str
cpu or gpu
- type: type
data type for calculations, default is np.float32
- alpha: array # TODO not implemented
frustration matrix
- num_sim: int
number of simulations
- method: str
integration method, default is heun
- same_initial_state: bool
use the same initial state for all simulations, default is False
- set_initial_state(nn, ns=1, engine='cpu', seed=None, same_initial_state=False, dtype=<class 'float'>)[source]¶
set initial state
- Parameters:
- nn: int
number of nodes
- ns: int
number of states
- engine: str
cpu or gpu
- seed: int
set random seed if not None
- same_initial_state: bool
use the same initial state for all simulations
- Returns:
- x: array [nn, ns]
initial state
vbi.models.cupy.mpr_modified_bold¶
- class MPR_sde(par: dict = {}, Bpar: dict = {})[source]¶
Montbrio-Pazo-Roxin model Cupy and Numpy implementation.
- Parameters:
- G: float. np.ndarray
global coupling strength
- dt: float
time step
- dt_bold: float
time step for Balloon model
- J: float, np.ndarray
model parameter
- eta: float, np.ndarray
model parameter
- tau:
model parameter
- delta:
model parameter
- tr: float
repetition time of fMRI
- noise_amp: float, np.array
amplitude of noise
- same_noise_per_sim:
same noise for all simulations
- iapp: float, np.ndarray
external input
- t_start: float
initial time
- t_cut: float
transition time
- t_end: float
end time
- num_nodes: int
number of nodes
- weights: np.ndarray
weighted connection matrix
- rv_decimate: int
sampling step from r and v
- output: str
output directory
- RECORD_TS: bool
store r and v time series
- num_sim: int
number of simulations
- method: str
integration method
- engine: str
cpu or gpu
- seed: int
seed for random number generator
- dtype: str
float or f
- initial_state: np.ndarray
initial state
- same_initial_state: bool
same initial state for all simulations
vbi.models.cupy.wilson_cowan¶
- class WC_sde(par: dict = {})[source]¶
Wilson-Cowan model of neural population dynamics.
Used Eqns 11 and 12 from [WC_1972] in the model dynamics. P and Q represent external inputs, which when exploring the phase portrait of the local model are set to constant values. However in the case of a full network, P and Q are the entry point to our long range and local couplings, that is, the activity from all other nodes is the external input to the local population [WC_1973], [D_2011].
References:
[WC_1972]Wilson, H.R. and Cowan, J.D. Excitatory and inhibitory interactions in localized populations of model neurons, Biophysical journal, 12: 1-24, 1972.
[WC_1973]Wilson, H.R. and Cowan, J.D A Mathematical Theory of the Functional Dynamics of Cortical and Thalamic Nervous Tissue
[D_2011]Daffertshofer, A. and van Wijk, B. On the influence of amplitude on the connectivity between phases Frontiers in Neuroinformatics, July, 2011
- prepare_input()[source]¶
Prepare input parameters, check dimensions and convert to cupy array if needed. Some parameters can be scalars, vectors or 2D arrays. 2D arrays parameters are heterogeneous accross nodes and simulations, but 1D arrays parameters are homogeneous accross nodes and heterogeneous accross simulations.
vector parameters: ns scalar parameters: 1 matrix parameters: nn x ns
vbi.models.cupy.ww¶
- class WW_sde(par: Dict = {}, Bpar: Dict = {})[source]¶
Wong-Wang neural mass including Excitatory and Inhibitory populations.
This model explicitly simulates both excitatory and inhibitory neural populations with distinct synaptic gating variables, firing rate dynamics, and network coupling. The model is based on the original Wong-Wang decision-making framework and has been extended for whole-brain network simulations.
- Main reference:
[original] Wong, K. F., & Wang, X. J. (2006). A recurrent network mechanism of time integration in perceptual decisions. Journal of Neuroscience, 26(4), 1314-1328.
- Additional references:
[reduced] Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann, P., & Corbetta, M. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. Journal of Neuroscience, 33(27), 11239-11252.
[original] Deco, G., Ponce-Alvarez, A., Hagmann, P., Romani, G. L., Mantini, D., & Corbetta, M. (2014). How local excitation-inhibition ratio impacts the whole brain dynamics. Journal of Neuroscience, 34(23), 7886-7898.
Parameters¶ Name
Explanation
Default Value
a_excExcitatory population gain parameter for firing rate transfer function (n/C).
310
a_inhInhibitory population gain parameter for firing rate transfer function (nC^-1).
0.615
b_excExcitatory population threshold parameter for firing rate transfer function (Hz).
125
b_inhInhibitory population threshold parameter for firing rate transfer function (Hz).
177
d_excExcitatory population saturation parameter for firing rate transfer function (s).
0.16
d_inhInhibitory population saturation parameter for firing rate transfer function (s).
0.087
tau_excExcitatory synaptic time constant (ms).
100.0
tau_inhInhibitory synaptic time constant (ms).
10.0
gamma_excExcitatory kinetic parameter (ms^-1).
0.641/1000.0
gamma_inhInhibitory kinetic parameter (ms^-1).
1.0/1000.0
W_excExcitatory population local weight for external input.
1.0
W_inhInhibitory population local weight for external input.
0.7
ext_currentExternal current input (nA). If array-like, it should be of shape (
num_nodes,num_sim).0.382
J_NMDANMDA synaptic coupling strength (nA).
0.15
J_IInhibitory synaptic coupling strength (nA).
1.0
w_plusLocal excitatory recurrence strength.
1.4
lambda_inh_excLong-range feedforward inhibition switch (0=off, 1=on).
0.0
G_excGlobal excitatory coupling strength. If array-like, it should be of length
num_sim.0.0
G_inhGlobal inhibitory coupling strength. If array-like, it should be of length
num_sim.0.0
t_endEnd time of simulation (ms).
1000.0
t_cutTime to cut off initial transient (ms).
0.0
dtTime step for integration (ms).
0.1
trRepetition time for BOLD signal (ms).
300.0
s_decimateDecimation factor for recording gating variables S.
1
sigmaNoise amplitude for stochastic integration.
0.0
weightsStructural connectivity matrix of shape (
num_nodes,num_nodes).None
num_simNumber of parallel simulations.
1
nnNumber of brain regions/nodes.
1
engineComputation engine: “cpu” or “gpu”.
“cpu”
dtypeData type: “float32” or “float”.
“float32”
seedRandom seed for reproducibility.
None
outputOutput directory for results.
“output”
initial_stateInitial state of the system of shape (2*`num_nodes`,
num_sim).None
same_initial_stateIf True, all simulations have the same initial state.
False
same_noise_per_simIf True, all simulations have the same noise realization.
False
RECORD_SIf True, record synaptic gating variables S.
False
RECORD_BOLDIf True, record BOLD signal using Balloon-Windkessel model.
True
- get_firing_rate(current: float, is_exc: bool = True)[source]¶
Calculate firing rate based on input current
- set_initial_state(nn, ns, engine, seed=None, same_initial_state=False, dtype=<class 'float'>)[source]¶
Set initial state
- Parameters:
- nnint
number of nodes
- nsint
number of simulations
- enginestr
cpu or gpu
- same_initial_conditionbool
same initial condition for all simulations
- seedint
random seed
- dtypestr
float: float64 f : float32
Numba models¶
vbi.models.numba.damp_oscillator¶
Damped Oscillator Model - Numba Implementation
This module implements a nonlinear damped oscillator model using Numba for high-performance numerical integration. The model exhibits rich dynamical behavior including fixed points, limit cycles, and chaotic dynamics depending on parameter values.
- The system is described by:
dx/dt = x - xy - ax²
dy/dt = xy - y - by²
Classes¶
DO_nb : Main model class Param : Parameter container (Numba jitclass)
Functions¶
euler, heun, rk4 : Integration methods _f_sys : System dynamics function _integrate : Main integration routine
- class Param(*args, **kwargs)[source]¶
Parameter container for the damped oscillator model.
A Numba-compiled parameter class that stores all model parameters for efficient access during numerical integration. This class is optimized for high-performance computing with Numba JIT compilation.
- Parameters:
- afloat, default 0.1
Damping parameter for x equation (ax² term).
- bfloat, default 0.05
Damping parameter for y equation (by² term).
- dtfloat, default 0.01
Integration time step.
- t_startfloat, default 0
Start time for integration.
- t_endfloat, default 100.0
End time for integration.
- t_cutfloat, default 20
Initial time to discard (burn-in period).
- outputstr, default “output”
Output directory name.
- methodstr, default “euler”
Integration method (“euler”, “heun”, or “rk4”).
- initial_statenp.ndarray, default [0.5, 1.0]
Initial conditions [x0, y0].
Notes
This class is compiled by Numba for efficient parameter access. The damped oscillator equations are:
\[ \begin{align}\begin{aligned}\frac{dx}{dt} = x - xy - ax^2\\\frac{dy}{dt} = xy - y - by^2\end{aligned}\end{align} \]
- euler(x, P)[source]¶
Euler integration method for the damped oscillator.
First-order explicit integration scheme: x(t+dt) = x(t) + dt * f(x, t)
- Parameters:
- xnp.ndarray
Current state vector [x, y].
- PParam
Parameter object containing dt and other parameters.
- Returns:
- np.ndarray
Next state vector after one integration step.
- heun(x, P)[source]¶
Heun’s integration method for the damped oscillator.
Second-order predictor-corrector method that provides better accuracy than Euler method with moderate computational cost.
- Parameters:
- xnp.ndarray
Current state vector [x, y].
- PParam
Parameter object containing dt and other parameters.
- Returns:
- np.ndarray
Next state vector after one integration step.
- rk4(x, P)[source]¶
Fourth-order Runge-Kutta integration method for the damped oscillator.
High-accuracy integration scheme that evaluates the system function four times per step for excellent precision.
- Parameters:
- xnp.ndarray
Current state vector [x, y].
- PParam
Parameter object containing dt and other parameters.
- Returns:
- np.ndarray
Next state vector after one integration step.
- class DO(par={})[source]¶
Numba implementation of the damped oscillator model.
The damped oscillator is a nonlinear dynamical system that exhibits complex behavior including limit cycles and chaotic dynamics. The system is described by the equations:
\[ \begin{align}\begin{aligned}\frac{dx}{dt} = x - xy - ax^2 \\\frac{dy}{dt} = xy - y - by^2\end{aligned}\end{align} \]where x and y are the state variables, and a and b are damping parameters.
Parameters¶ Name
Explanation
Default Value
aDamping parameter for the first state variable x.
0.1
bDamping parameter for the second state variable y.
0.05
dtIntegration time step.
0.01
t_startStart time of simulation.
0
t_endEnd time of simulation.
100.0
t_cutTime from which to start collecting output (burn-in period).
20
outputOutput directory path for saving results.
“output”
methodIntegration method. Options: “euler”, “heun”, “rk4”.
“euler”
initial_stateInitial state vector [x0, y0] of length 2.
[0.5, 1.0]
- Usage example:
Notes
The damped oscillator model is a classic example of a nonlinear dynamical system. Depending on the parameter values, the system can exhibit:
Fixed points (stable equilibria)
Limit cycles (periodic oscillations)
Chaotic behavior (sensitive dependence on initial conditions)
The model supports three integration methods:
Euler: First-order explicit method, fastest but least accurate
Heun: Second-order predictor-corrector method, good balance of speed and accuracy
RK4: Fourth-order Runge-Kutta method, most accurate but slowest
The simulation includes a burn-in period (
t_cut) to allow transients to decay before collecting output data.
- __init__(par={})[source]¶
Initialize the damped oscillator model.
- Parameters:
- pardict, optional
Dictionary containing model parameters. See class documentation for available parameters. Default is an empty dictionary which uses all default parameter values.
- check_parameters(par)[source]¶
Validate model parameters.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter name is not recognized.
- get_parobj(par={})[source]¶
Create parameter object with default values and user overrides.
- Parameters:
- pardict, optional
Dictionary of parameters to override defaults.
- Returns:
- Param
Numba jitclass parameter object.
- update_par(par={})[source]¶
Update model parameters.
- Parameters:
- pardict, optional
Dictionary of parameters to update. Keys must be valid parameter names.
- run(par={}, x0=None, verbose=False)[source]¶
Run the damped oscillator simulation.
- Parameters:
- pardict, optional
Dictionary of parameters to update for this simulation run. Any parameter from the class documentation can be updated.
- x0array-like, optional
Initial state vector [x0, y0] of length 2. If None, uses the initial state set during initialization.
- verbosebool, optional
If True, print verbose output during simulation. Default is False.
- Returns:
- dict
Dictionary containing simulation results:
‘t’ : np.ndarray of shape (n_steps,) - time points
‘x’ : np.ndarray of shape (n_steps, 2) - simulated trajectories where x[:, 0] is x(t) and x[:, 1] is y(t)
vbi.models.numba.ghb¶
- run(P, times)[source]¶
Run the Generic Hopf Bifurcation model simulation.
This function integrates the coupled Generic Hopf oscillators with BOLD signal generation using Euler-Maruyama stochastic integration.
- Parameters:
- PParGHB
Parameter object containing all model parameters.
- timesnp.ndarray
Time array for the simulation.
- Returns:
- tuple
(t_bold, bold) where t_bold is time array and bold is BOLD signal of shape (nn, n_timepoints).
- class GHB_sde(par: dict = {})[source]¶
Numba implementation of the Generic Hopf Bifurcation model with BOLD signal generation.
This model implements a network of coupled Stuart-Landau oscillators (Generic Hopf Bifurcation) with additive noise and hemodynamic response modeling to generate BOLD signals. Each brain region is modeled as a 2D oscillator with intrinsic frequency and bifurcation parameter.
The neural dynamics follow: dx/dt = x(η - x² - y²) - ωy + G·coupling_x + noise dy/dt = y(η - x² - y²) + ωx + G·coupling_y + noise
Where: - (x, y) are the neural state variables in Cartesian coordinates - η (eta) is the bifurcation parameter controlling oscillation amplitude - ω (omega) is the intrinsic frequency of oscillation - G is the global coupling strength - The (x² + y²) term provides amplitude-dependent damping
Parameters¶ Name
Explanation
Default Value
GGlobal coupling strength scaling network connections.
1.0
dtIntegration time step in seconds.
0.001
sigmaNoise amplitude (standard deviation of additive Gaussian noise).
0.1
tendEnd time of simulation in seconds.
10.0
tcutTime from which to start collecting output (burn-in period).
0.0
etaBifurcation parameter array of length
nn. Controls oscillation amplitude. If array-like, it should be of lengthnn(number of nodes).np.array([])
omegaIntrinsic frequency array of length
nn. Sets the natural oscillation frequency for each region. If array-like, it should be of lengthnn.np.array([])
weightsStructural connectivity matrix of shape (
nn,nn). Must be provided.np.array([[], []])
decimateDecimation factor for output time series (every
decimate-th point is saved).1
init_stateInitial state vector of shape (2*nn,) containing [x₀, y₀]. If None, random initial conditions are generated.
np.array([])
seedRandom seed for reproducible simulations. If -1, no seeding is applied.
-1
- Usage example:
Notes
The model combines neural oscillator dynamics with a simplified BOLD hemodynamic response based on the Balloon-Windkessel model. The BOLD signal is computed using the Friston 2003 formulation with fixed parameters optimized for realistic hemodynamic responses.
The Generic Hopf Bifurcation model is particularly useful for studying: - Oscillatory brain dynamics and synchronization - Transitions between different dynamical regimes - Frequency-specific network interactions - BOLD signal generation from neural oscillations
- __init__(par: dict = {}) None[source]¶
Initialize the Generic Hopf Bifurcation model.
- Parameters:
- pardict, optional
Dictionary containing model parameters. See class documentation for available parameters. The ‘weights’, ‘eta’, and ‘omega’ parameters are required for proper functionality.
- get_par_obj(par: dict)[source]¶
Create parameter object from dictionary.
- Parameters:
- pardict
Parameter dictionary.
- Returns:
- ParGHB
Numba jitclass parameter object.
- check_parameters(par: dict) None[source]¶
Validate that all provided parameters are recognized.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter name is not recognized.
- set_initial_state(seed=None)[source]¶
Generate random initial conditions for the model.
- Parameters:
- seedint, optional
Random seed for reproducible initial conditions.
- Returns:
- np.ndarray
Initial state vector of shape (2*nn,) with random values in [0, 1].
- check_input()[source]¶
Validate model parameters before simulation.
- Raises:
- AssertionError
If weights matrix is None, not square, or if eta/omega arrays don’t match the number of nodes.
- run(par={}, tspan=None, x0=None, verbose=True)[source]¶
Run the Generic Hopf Bifurcation simulation.
- Parameters:
- pardict, optional
Dictionary of parameters to update for this simulation run. Any parameter from the class documentation can be updated.
- tspantuple or None, optional
Time span as (start_time, end_time). If None, uses (0, tend).
- x0np.ndarray, optional
Initial state vector of shape (2*nn,). If None, random initial conditions are generated using set_initial_state().
- verbosebool, optional
Whether to print verbose output (currently unused).
- Returns:
- dict
Dictionary containing simulation results:
‘t’: np.ndarray of time points (after burn-in and decimation)
‘bold’: np.ndarray of shape (nn, n_timepoints) containing simulated BOLD signals for each brain region
- class ParGHB(*args, **kwargs)[source]¶
Numba jitclass container for Generic Hopf Bifurcation model parameters.
This class holds all parameters needed for the Generic Hopf Bifurcation model in a format optimized for Numba compilation. It stores both scalar parameters and array parameters like connectivity weights, bifurcation parameters, and frequencies.
Note: This is an internal class used by the GHB_sde class. Users should not instantiate this class directly.
vbi.models.numba.jansen_rit¶
- set_initial_state_jr(nn, seed=-1)[source]¶
Initial state for JR: stack 6*n vectors. Mirrors ranges similar to the CuPy reference implementation.
- class ParJR(*args, **kwargs)[source]¶
Numba jitclass container for Jansen-Rit model parameters.
This class holds all parameters needed for the Jansen-Rit neural mass model in a format optimized for Numba compilation. It stores both scalar parameters and array parameters like connectivity weights and coupling constants.
Note: This is an internal class used by the JR_sde class. Users should not instantiate this class directly.
- f_jr(x, t, P)[source]¶
Compute the right-hand side of the Jansen-Rit differential equations.
This function implements the deterministic part of the Jansen-Rit neural mass model equations for a network of coupled brain regions.
- Parameters:
- xnp.ndarray
Current state vector of shape (6*nn,) containing stacked arrays: [x, y, z, x’, y’, z’] where nn is the number of nodes.
- tfloat
Current time (not used in autonomous system).
- PParJR
Parameter object containing all model parameters.
- Returns:
- np.ndarray
Derivative vector of shape (6*nn,) containing dx/dt.
- heun_sde(x, t, P)[source]¶
Perform one step of Heun’s method for stochastic differential equations.
This implements the Heun scheme (also called improved Euler method) for integrating stochastic differential equations with additive noise.
- Parameters:
- xnp.ndarray
Current state vector of shape (6*nn,).
- tfloat
Current time.
- PParJR
Parameter object containing dt, sigma, and other model parameters.
- Returns:
- np.ndarray
Updated state vector after one integration step.
- integrate(P)[source]¶
Integrate the Jansen-Rit model over time.
This function performs the main time integration loop for the Jansen-Rit model, using the Heun stochastic integration scheme. It includes options for burn-in period and decimation of the output.
- Parameters:
- PParJR
Parameter object containing all simulation parameters including initial_state, t_end, dt, t_cut, and decimate.
- Returns:
- dict
Dictionary with keys:
‘t’: np.ndarray of time points (decimated)
‘x’: np.ndarray of shape (n_steps, nn) containing the output time series (y - z) representing local field potentials
- class JR_sde(par_jr: dict)[source]¶
Numba implementation of the Jansen-Rit neural mass model.
Parameters¶ Name
Explanation
Default Value
AExcitatory post synaptic potential amplitude.
3.25
BInhibitory post synaptic potential amplitude.
22.0
aInverse time constant of the excitatory postsynaptic potential (1/a = time constant).
0.1 (time constant: 10.0)
bInverse time constant of the inhibitory postsynaptic potential (1/b = time constant).
0.05 (time constant: 20.0)
C0Average number of synapses between pyramidal cells and excitatory interneurons. If array-like, it should be of length
nn(number of nodes).135.0
C1Average number of synapses between excitatory interneurons and pyramidal cells. If array-like, it should be of length
nn.0.8 * 135.0
C2Average number of synapses between pyramidal cells and inhibitory interneurons. If array-like, it should be of length
nn.0.25 * 135.0
C3Average number of synapses between inhibitory interneurons and pyramidal cells. If array-like, it should be of length
nn.0.25 * 135.0
vmaxMaximum firing rate of the sigmoid function.
0.005
v0Potential at half of maximum firing rate (inflection point of sigmoid).
6.0
rSlope of sigmoid function at
v0.0.56
GGlobal coupling strength scaling the network connections.
1.0
muMean input to the excitatory population (external drive).
0.24
noise_ampAmplitude of the stochastic noise applied to the excitatory population.
0.01
weightsStructural connectivity matrix of shape (
nn,nn). Must be provided.None
dtIntegration time step.
0.01
t_endEnd time of simulation.
1000.0
t_cutTime from which to start collecting output (burn-in period).
0.0
decimateDecimation factor for the output time series (every
decimate-th point is saved).1
seedRandom seed for reproducible simulations. If -1 or None, no seeding is applied.
-1
initial_stateInitial state vector of shape (6*nn,). If None, random initial conditions are generated.
None
- Usage example (single simulation):
Notes
The Jansen-Rit model describes the dynamics of a cortical column with three neural populations: - Pyramidal cells (main excitatory population) - Excitatory interneurons - Inhibitory interneurons
The model equations are integrated using the Heun stochastic integration scheme. The output represents the difference between excitatory and inhibitory postsynaptic potentials (y - z), which corresponds to the local field potential that can be measured experimentally.
- __init__(par_jr: dict)[source]¶
Initialize the Jansen-Rit model.
- Parameters:
- par_jrdict
Dictionary containing model parameters. See class documentation for available parameters. The ‘weights’ parameter is required and must be a square connectivity matrix.
- check_parameters(par: dict)[source]¶
Validate provided parameter names.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter names are invalid.
- check_input()[source]¶
Validate model parameters.
- Raises:
- ValueError
If any parameter values are invalid (e.g., t_cut >= t_end, decimate < 1, or dimension mismatches).
- set_initial_state(seed: int | None = None)[source]¶
Set random initial state for the simulation.
- Parameters:
- seedint, optional
Random seed for reproducible initial conditions. If None, uses the seed specified during initialization.
- run(par: dict | None = None, x0: ndarray | None = None)[source]¶
Run the Jansen-Rit simulation.
- Parameters:
- pardict, optional
Dictionary of parameters to update for this simulation run. Any parameter from the class documentation can be updated.
- x0np.ndarray, optional
Initial state vector of shape (6*nn,). If None, uses the initial state set during initialization or by set_initial_state().
- Returns:
- dict
Dictionary containing simulation results:
‘t’: np.ndarray of shape (n_steps,) - time points
‘x’: np.ndarray of shape (n_steps, nn) - simulated time series (y - z) representing local field potentials
vbi.models.numba.mpr¶
- f_mpr(x, t, P)[source]¶
Compute the right-hand side of the Montbrió neural mass model.
This function implements the deterministic part of the Montbrió model equations for a network of coupled neural populations. The model describes the macroscopic dynamics of quadratic integrate-and-fire (QIF) neuron populations in terms of population firing rate (r) and mean membrane potential (v).
The equations are: τ dr/dt = Δ/(τπ) + 2rv τ dv/dt = v² + η + I_app + Jr - (πτr)² + G·coupling
- Parameters:
- xnp.ndarray
Current state vector of shape (2*nn,) containing stacked arrays: [r₀, r₁, …, rₙ, v₀, v₁, …, vₙ] where nn is the number of nodes. r represents the population firing rate and v represents the mean membrane potential.
- tfloat
Current time (not used in autonomous system).
- PParMPR
Parameter object containing all model parameters.
- Returns:
- np.ndarray
Derivative vector of shape (2*nn,) containing dx/dt.
- heun_sde(x, t, P)[source]¶
Perform one step of Heun’s method for stochastic differential equations.
This implements the Heun scheme for integrating the Montbrió model with additive noise applied to both firing rate and membrane potential variables.
- Parameters:
- xnp.ndarray
Current state vector of shape (2*nn,).
- tfloat
Current time.
- PParMPR
Parameter object containing dt, sigma_r, sigma_v and other parameters.
- Returns:
- np.ndarray
Updated state vector after one integration step.
- do_bold_step(r_in, s, f, ftilde, vtilde, qtilde, v, q, dtt, P)[source]¶
Perform one step of BOLD signal computation using the Balloon-Windkessel model.
This function implements the hemodynamic response model that converts neural activity (firing rate) into BOLD signal through a cascade of physiological processes including blood flow, volume, and oxygenation changes.
- Parameters:
- r_infloat
Input neural activity (firing rate).
- s, f, ftilde, vtilde, qtilde, v, qnp.ndarray
BOLD state variables (flow, volume, deoxygenation).
- dttfloat
Time step for BOLD integration.
- PParBold
BOLD parameter object.
- integrate(P, B)[source]¶
Integrate the Montbrió model over time with BOLD signal computation.
This function performs the main time integration loop for the Montbrió model, using the Heun stochastic integration scheme. It optionally computes both neural activity (r, v) and hemodynamic BOLD signals.
- Parameters:
- PParMPR
Parameter object containing all simulation parameters.
- BParBold
BOLD parameter object for hemodynamic response computation.
- Returns:
- dict
Dictionary with keys: - ‘rv_t’: np.ndarray of time points for neural activity - ‘rv_d’: np.ndarray of shape (n_steps, 2*nn) containing [r, v] time series - ‘bold_t’: np.ndarray of time points for BOLD signal - ‘bold_d’: np.ndarray of shape (n_steps, nn) containing BOLD time series
- class MPR_sde(par_mpr: dict = {})[source]¶
Numba implementation of the Montbrió neural mass model with BOLD signal generation.
This model implements the exact macroscopic dynamics derived from infinitely all-to-all coupled quadratic integrate-and-fire (QIF) neurons in the thermodynamic limit. The model describes the collective firing activity (r) and mean membrane potential (v) of neural populations, providing a rigorous mathematical foundation for whole-brain network modeling.
The neural dynamics follow the Montbrió equations: τ dr/dt = 2rv + Δ/(πτ) τ dv/dt = v² - (πτr)² + Jτr + η + G·coupling + I_stim + noise
Where: - r is the population firing rate - v is the mean membrane potential - τ is the characteristic time constant - J is the synaptic weight - Δ is the spread of heterogeneous excitabilities - η is the mean excitability - G is the global coupling strength
The model can operate in bistable regime, exhibiting down-state (low firing) and up-state (high firing) attractors, which is fundamental for realistic brain dynamics and functional connectivity patterns.
Parameters¶ Name
Explanation
Default Value
GGlobal coupling strength scaling network connections.
0.5
dtIntegration time step in milliseconds.
0.01
JSynaptic weight strength in ms⁻¹.
14.5
etaMean excitability parameter. If array-like, it should be of length
nn.np.array([-4.6])
tauCharacteristic time constant in ms.
1.0
deltaSpread of heterogeneous excitability distribution in ms⁻¹.
0.7
rv_decimateDecimation factor for neural activity recording.
1.0
noise_ampAmplitude of additive Gaussian noise.
0.037
weightsStructural connectivity matrix of shape (
nn,nn). Must be provided.np.array([[], []])
t_initInitial time of simulation.
0.0
t_cutTime from which to start collecting output (burn-in period).
0.0
t_endEnd time of simulation in milliseconds.
1000.0
iappExternal applied current (stimulation).
0.0
seedRandom seed for reproducible simulations. If -1, no seeding is applied.
-1
outputOutput directory name.
“output”
RECORD_RVWhether to record neural activity (r, v) time series.
True
RECORD_BOLDWhether to record BOLD signal time series.
True
trBOLD repetition time in milliseconds.
500.0
- Usage example:
Notes
The Montbrió model provides an exact mathematical description derived from microscopic spiking neuron networks, making it particularly suitable for: - Whole-brain network modeling with rigorous theoretical foundation - Studying bistable dynamics and metastable states - Investigating the relationship between neural activity and BOLD signals - Parameter estimation and model validation against empirical data
The model includes sophisticated BOLD signal generation using the Balloon-Windkessel hemodynamic response model, enabling direct comparison with fMRI measurements.
References
Montbrió, E., et al. (2015). Macroscopic description for networks of spiking neurons. Physical Review X, 5(2), 021028.
- __init__(par_mpr: dict = {}) None[source]¶
Initialize the Montbrió model.
- Parameters:
- par_mprdict, optional
Dictionary containing model parameters. See class documentation for available parameters. The ‘weights’ parameter is required for proper functionality.
- check_parameters(par: dict) None[source]¶
Validate that all provided parameters are recognized.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter name is not recognized.
- get_par_mpr(par: dict)[source]¶
Create parameter object from dictionary with validation.
- Parameters:
- pardict
Parameter dictionary.
- Returns:
- ParMPR
Numba jitclass parameter object with validated parameters.
- set_initial_state()[source]¶
Generate random initial conditions for the model.
Sets the initial state for both firing rate (r) and membrane potential (v) variables using random values appropriate for the Montbrió model dynamics.
- check_input()[source]¶
Validate model parameters before simulation.
- Raises:
- AssertionError
If weights matrix is None, not square, if initial_state length doesn’t match 2*nn, or if t_cut >= t_end.
- run(par={}, x0=None, verbose=True)[source]¶
Run the Montbrió simulation.
- Parameters:
- pardict, optional
Dictionary of parameters to update for this simulation run. Any parameter from the class documentation can be updated.
- x0np.ndarray, optional
Initial state vector of shape (2*nn,) containing [r₀, v₀]. If None, random initial conditions are generated.
- verbosebool, optional
Whether to print verbose output (currently unused).
- Returns:
- dict
Dictionary containing simulation results:
‘rv_t’: np.ndarray of time points for neural activity (in ms)
‘rv_d’: np.ndarray of shape (n_steps, 2*nn) containing [r, v] time series (firing rate and membrane potential)
‘bold_t’: np.ndarray of time points for BOLD signal (in ms)
‘bold_d’: np.ndarray of shape (n_steps, nn) containing simulated BOLD signals for each brain region
- set_initial_state(nn, seed=None)[source]¶
Generate random initial state for the Montbrió model.
Creates initial conditions with firing rates in [0, 1.5] and membrane potentials in [-2, 2], appropriate for the model dynamics.
- Parameters:
- nnint
Number of nodes/brain regions.
- seedint, optional
Random seed for reproducible initial conditions.
- Returns:
- np.ndarray
Initial state vector of shape (2*nn,) with [r₀, v₀] values.
- class ParMPR(*args, **kwargs)[source]¶
Numba jitclass container for Montbrió model parameters.
This class holds all parameters needed for the Montbrió neural mass model in a format optimized for Numba compilation. It stores both scalar parameters and array parameters like connectivity weights and excitability values.
Note: This is an internal class used by the MPR_sde class. Users should not instantiate this class directly.
- class ParBold(*args, **kwargs)[source]¶
Numba jitclass container for BOLD hemodynamic model parameters.
This class holds parameters for the Balloon-Windkessel model used to convert neural activity into simulated BOLD signals. Based on the Friston 2003 hemodynamic response model.
Note: This is an internal class used by the MPR_sde class. Users should not instantiate this class directly.
vbi.models.numba.vep¶
- seed_rng(seed: int)[source]¶
Seed the NumPy random number generator for reproducible results.
- Parameters:
- seedint
Random seed value. If negative, no seeding is performed.
- class ParVEP(*args, **kwargs)[source]¶
Parameter container class for VEP model (Numba jitclass).
This class holds all parameters needed for the VEP simulation and is compiled by Numba for efficient access during integration. It automatically calculates the number of nodes (nn) from the weights matrix shape.
- Parameters:
- Gfloat, optional
Global coupling strength (default: 1.0)
- dtfloat, optional
Integration time step in seconds (default: 0.01)
- taufloat, optional
Time constant for slow variable dynamics (default: 10.0)
- etaarray_like, optional
Epileptogenicity per region (default: [-1.5])
- iextarray_like, optional
External input per region (default: [0.0])
- weightsarray_like
Connectivity matrix (nn x nn)
- t_cutfloat, optional
Time to cut initial transient (default: 0.0)
- t_endfloat, optional
End time of simulation (default: 100.0)
- methodstr, optional
Integration method: “euler” or “heun” (default: “euler”)
- seedint, optional
Random seed for noise generation (default: -1)
- initial_statearray_like, optional
Initial state vector (default: [0.0, 0.0])
- sigmafloat, optional
Noise amplitude (default: 0.1)
- record_stepint, optional
Recording decimation factor (default: 1)
- outputstr, optional
Output directory (default: “output”)
- f_vep(x, t, P)[source]¶
Right-hand side of VEP dynamics equations.
Computes the time derivatives for the VEP model: - dx/dt = 1 - x³ - 2x² - y + I_ext - dy/dt = (4(x - η) - y - G*coupling) / τ
- Parameters:
- xndarray
State vector of shape (2*nn,) where: - x[0:nn] contains x-variables (fast variables) - x[nn:2*nn] contains y-variables (slow variables)
- tfloat
Current time (unused but required for integrator interface)
- PParVEP
Parameter object containing model parameters
- Returns:
- ndarray
Time derivatives of shape (2*nn,)
- euler_step(x, t, P)[source]¶
Perform one Euler integration step with additive noise.
- Parameters:
- xndarray
Current state vector of shape (2*nn,)
- tfloat
Current time
- PParVEP
Parameter object containing model parameters
- Returns:
- ndarray
Updated state vector after one integration step
- heun_step(x, t, P)[source]¶
Perform one Heun integration step with additive noise.
Heun’s method is a second-order Runge-Kutta method that provides better accuracy than Euler’s method for stochastic differential equations.
- Parameters:
- xndarray
Current state vector of shape (2*nn,)
- tfloat
Current time
- PParVEP
Parameter object containing model parameters
- Returns:
- ndarray
Updated state vector after one integration step
- set_initial_state_jit(nn: int, seed: int)[source]¶
Generate random initial state for VEP model.
Creates initial conditions matching the C++/Python implementation: - Fast variables (x): uniformly distributed in [-3, -2] - Slow variables (y): uniformly distributed in [0, 3.5]
- Parameters:
- nnint
Number of brain regions/nodes
- seedint
Random seed for reproducibility. If negative, no seeding is performed.
- Returns:
- ndarray
Initial state vector of shape (2*nn,) with: - x0[:nn] ~ U(-3, -2) (fast variables) - x0[nn:] ~ U(0, 3.5) (slow variables)
- class VEP_sde(par_vep: dict = {})[source]¶
Virtual Epileptic Patient (VEP) model - Numba implementation.
The VEP model is a 2D reduction of the full Epileptor model, designed for personalized whole-brain network modeling of epilepsy spread. This model provides a comprehensive description of epileptic seizures and has been extensively used in clinical applications for seizure prediction and understanding epilepsy dynamics.
The model equations are:
\[\begin{split}\frac{dx_i}{dt} &= 1 - x_i^3 - 2x_i^2 - y_i + I_{ext,i} \\ \frac{dy_i}{dt} &= \frac{1}{\tau}(4(x_i - \eta_i) - y_i - G \sum_j W_{ij}(x_j - x_i))\end{split}\]where \(x_i\) and \(y_i\) are the fast and slow variables at region \(i\), \(\eta_i\) represents the epileptogenicity, and the network coupling uses a Laplacian form with connectivity matrix \(W_{ij}\).
- Main references:
Jirsa, V.K., et al. (2017). The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread. NeuroImage, 145, 377-388.
Jirsa, V.K., et al. (2014). On the nature of seizure dynamics. Brain, 137(8), 2210-2230.
Parameters¶ Name
Explanation
Default Value
GGlobal coupling strength that scales network interactions.
1.0
dtTime step for numerical integration (s).
0.01
tauTime constant for the slow variable dynamics (s).
10.0
etaEpileptogenicity parameter per region. If scalar, broadcasted to all regions.
-1.5
iextExternal input current per region (nA). If scalar, broadcasted to all regions.
0.0
weightsStructural connectivity matrix of shape (
nn,nn). Required parameter.None
t_cutTime to discard initial transient (s).
0.0
t_endTotal simulation time (s).
100.0
methodIntegration method: “euler” or “heun”.
“euler”
seedRandom seed for initial state generation. Use -1 for random initialization.
-1
initial_stateInitial state vector of shape (2*`nn`,). If None, random state is generated.
None
sigmaNoise amplitude for stochastic integration.
0.1
record_stepDecimation factor for recording output (saves every nth step).
1
outputOutput directory path (currently unused in this implementation).
“output”
- Returns:
- dict
Dictionary with keys: - ‘t’: 1D array of time points (after t_cut, with optional decimation) - ‘x’: 2D array of shape (
nn,nt_saved) containing only the fast variableNotes
The initial state is automatically generated as: x[:nn] ~ U(-3, -2), y[nn:] ~ U(0, 3.5)
Network coupling uses Laplacian form: \(\sum_j W_{ij}(x_j - x_i)\)
Noise is additive with strength \(\sigma \sqrt{dt}\) applied to all state variables
Only the fast variable (x) is recorded in the output, matching clinical data processing
Examples
>>> import numpy as np >>> from vbi.models.numba.vep import VEP_sde >>> nn = 4 >>> W = (np.ones((nn, nn)) - np.eye(nn)) * 0.5 >>> params = { ... "G": 1.0, ... "seed": 42, ... "weights": W, ... "tau": 10.0, ... "eta": -3.5, ... "sigma": 0.0, ... "iext": 3.1, ... "dt": 0.1, ... "t_end": 14.0, ... "t_cut": 1.0, ... "record_step": 1, ... "method": "heun", ... } >>> model = VEP_sde(params) >>> result = model.run() >>> t = result['t'] >>> x = result['x'] >>> print(f"Time shape: {t.shape}, Data shape: {x.shape}") Time shape: (130,), Data shape: (4, 130)
vbi.models.numba.wilson_cowan¶
- check_vec_size(x, nn)[source]¶
Return a length-nn vector from scalar/1-vector or already-length-nn input.
- class ParWC(*args, **kwargs)[source]¶
Parameter class for the Wilson-Cowan neural mass model.
This Numba jitclass holds all parameters required for Wilson-Cowan simulation including synaptic strengths, time constants, sigmoid parameters, noise levels, and simulation settings. Uses Numba for high-performance compilation.
- Parameters:
- c_eearray-like, default [16.0]
Excitatory to excitatory synaptic strength
- c_eiarray-like, default [12.0]
Excitatory to inhibitory synaptic strength
- c_iearray-like, default [15.0]
Inhibitory to excitatory synaptic strength
- c_iiarray-like, default [3.0]
Inhibitory to inhibitory synaptic strength
- tau_earray-like, default [8.0]
Excitatory population time constant
- tau_iarray-like, default [8.0]
Inhibitory population time constant
- a_efloat, default 1.3
Excitatory sigmoid maximum firing rate parameter
- a_ifloat, default 2.0
Inhibitory sigmoid maximum firing rate parameter
- b_efloat, default 4.0
Excitatory sigmoid steepness parameter
- b_ifloat, default 3.7
Inhibitory sigmoid steepness parameter
- c_efloat, default 1.0
Excitatory sigmoid center parameter
- c_ifloat, default 1.0
Inhibitory sigmoid center parameter
- theta_efloat, default 0.0
Excitatory threshold parameter
- theta_ifloat, default 0.0
Inhibitory threshold parameter
- r_efloat, default 1.0
Excitatory refractory parameter
- r_ifloat, default 1.0
Inhibitory refractory parameter
- k_efloat, default 0.994
Excitatory maximum response parameter
- k_ifloat, default 0.999
Inhibitory maximum response parameter
- alpha_efloat, default 1.0
Excitatory scaling parameter
- alpha_ifloat, default 1.0
Inhibitory scaling parameter
- Parray-like, default [0.0]
External input to excitatory population
- Qarray-like, default [0.0]
External input to inhibitory population
- g_efloat, default 0.0
Excitatory global coupling strength
- g_ifloat, default 0.0
Inhibitory global coupling strength
- dtfloat, default 0.01
Integration time step (ms)
- t_endfloat, default 300.0
Simulation end time (ms)
- t_cutfloat, default 0.0
Initial time to discard (ms)
- sigmoid_vec(x, a, b, c, shift_sigmoid)[source]¶
Vectorized sigmoidal transfer function for Wilson-Cowan model.
Computes either standard sigmoid or shifted sigmoid that maps synaptic input to firing rate, capturing saturation effects and thresholds.
- Parameters:
- xnp.ndarray
Input values (synaptic drive).
- afloat
Sigmoid slope parameter.
- bfloat
Sigmoid threshold parameter.
- cfloat
Maximum output of sigmoid.
- shift_sigmoidbool
If True, uses shifted sigmoid: c * (sigmoid(a(x-b)) - sigmoid(-ab)) If False, uses standard sigmoid: c / (1 + exp(-a(x-b)))
- Returns:
- np.ndarray
Sigmoid-transformed firing rates.
- f_wc(x, t, P)[source]¶
Compute the right-hand side of the Wilson-Cowan neural mass model.
This function implements the deterministic part of the Wilson-Cowan equations for a network of coupled excitatory-inhibitory neural populations. The model describes the temporal evolution of mean firing rates using nonlinear differential equations with sigmoidal transfer functions.
The equations are:
\[ \begin{align}\begin{aligned}\tau_e \frac{dE}{dt} = -E + (k_e - r_e E) \cdot S_e(\alpha_e(c_{ee}E - c_{ei}I + P - \theta_e + g_e\sum w_{ij}E_j))\\\tau_i \frac{dI}{dt} = -I + (k_i - r_i I) \cdot S_i(\alpha_i(c_{ie}E - c_{ii}I + Q - \theta_i + g_i\sum w_{ij}I_j))\end{aligned}\end{align} \]
- Parameters:
- xnp.ndarray
Current state vector of shape (2*nn,) containing stacked arrays: [E₀, E₁, …, Eₙ, I₀, I₁, …, Iₙ] where nn is the number of nodes. E represents excitatory population activity, I represents inhibitory.
- tfloat
Current time (not used in autonomous system).
- PParWC
Parameter object containing all model parameters.
- Returns:
- np.ndarray
Derivative vector of shape (2*nn,) containing dx/dt.
- class WC_sde_numba(par: dict = {})[source]¶
Numba implementation of the Wilson-Cowan neural mass model with stochastic dynamics.
The Wilson-Cowan model is a seminal neural mass model that describes the dynamics of connected excitatory and inhibitory neural populations at the cortical microcolumn level. It models the temporal evolution of mean firing rates using nonlinear differential equations with sigmoidal transfer functions that capture saturation effects and thresholds in neural response.
This implementation includes: - Coupled excitatory (E) and inhibitory (I) populations per brain region - Network connectivity through structural connectivity matrix - Additive Gaussian noise for stochastic dynamics - Flexible sigmoid functions (standard or shifted) - Efficient Numba compilation for high-performance simulation
The model equations are:
\[ \begin{align}\begin{aligned}\tau_e \frac{dE}{dt} = -E + (k_e - r_e E) \cdot S_e(\alpha_e(c_{ee}E - c_{ei}I + P - \theta_e + g_e\sum w_{ij}E_j)) + \text{noise}\\\tau_i \frac{dI}{dt} = -I + (k_i - r_i I) \cdot S_i(\alpha_i(c_{ie}E - c_{ii}I + Q - \theta_i + g_i\sum w_{ij}I_j)) + \text{noise}\end{aligned}\end{align} \]
Parameters¶ Name
Explanation
Default Value
c_eeExcitatory to excitatory synaptic strength. If array-like, should be of length
nn.16.0
c_eiInhibitory to excitatory synaptic strength. If array-like, should be of length
nn.12.0
c_ieExcitatory to inhibitory synaptic strength. If array-like, should be of length
nn.15.0
c_iiInhibitory to inhibitory synaptic strength. If array-like, should be of length
nn.3.0
tau_eTime constant of excitatory population. If array-like, should be of length
nn.8.0
tau_iTime constant of inhibitory population. If array-like, should be of length
nn.8.0
a_eSigmoid slope for excitatory population.
1.3
a_iSigmoid slope for inhibitory population.
2.0
b_eSigmoid threshold for excitatory population.
4.0
b_iSigmoid threshold for inhibitory population.
3.7
c_eMaximum output of sigmoid for excitatory population.
1.0
c_iMaximum output of sigmoid for inhibitory population.
1.0
theta_eFiring threshold for excitatory population.
0.0
theta_iFiring threshold for inhibitory population.
0.0
r_eRefractoriness of excitatory population.
1.0
r_iRefractoriness of inhibitory population.
1.0
k_eScaling constant for excitatory output.
0.994
k_iScaling constant for inhibitory output.
0.999
alpha_eGain of excitatory population.
1.0
alpha_iGain of inhibitory population.
1.0
PExternal input to excitatory population. If array-like, should be of length
nn.0.0
QExternal input to inhibitory population. If array-like, should be of length
nn.0.0
g_eGlobal coupling strength for excitatory populations.
0.0
g_iGlobal coupling strength for inhibitory populations.
0.0
weightsStructural connectivity matrix of shape (
nn,nn). Must be provided.None
dtIntegration time step.
0.01
t_endEnd time of simulation.
300.0
t_cutTime from which to start collecting output (burn-in period).
0.0
noise_ampAmplitude of additive Gaussian noise.
0.0
decimateDecimation factor for output time series (every
decimate-th point is saved).1
RECORD_EIWhich populations to record: “E” (excitatory only), “I” (inhibitory only), “EI” (both).
“E”
shift_sigmoidWhether to use shifted sigmoid function.
False
seedRandom seed for reproducible simulations. If -1, no seeding is applied.
-1
initial_stateInitial state vector of shape (2*nn,). If None, random initial conditions are generated.
None
- Usage example:
Notes
The Wilson-Cowan model is widely used for understanding: - Oscillations and wave propagation in neural tissue - Pattern formation and spatial dynamics - Responses to external stimuli and perturbations - Brain dysfunction in neurological conditions (e.g., Parkinson’s disease) - Local field potentials (LFPs) and EEG signal generation
The model’s nonlinear dynamics arise from the sigmoidal transfer functions that map synaptic input to firing rate, allowing for rich dynamical behaviors including fixed points, limit cycles, and complex spatiotemporal patterns.
References
Wilson, H. R., & Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical journal, 12(1), 1-24.
- __init__(par: dict = {})[source]¶
Initialize the Wilson-Cowan model.
- Parameters:
- pardict, optional
Dictionary containing model parameters. See class documentation for available parameters. The ‘weights’ parameter is required.
- set_initial_state()[source]¶
Generate random initial conditions for the model.
Sets random initial state for both excitatory and inhibitory populations with values appropriate for Wilson-Cowan dynamics.
- check_parameters(par: dict) None[source]¶
Validate that all provided parameters are recognized.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter name is not recognized.
- check_input()[source]¶
Check shape and consistency of input parameters.
- Returns:
- bool
True if input is provided, False otherwise.
- run(par: dict | None = None, x0=None, verbose: bool = True)[source]¶
Run the Wilson-Cowan model simulation.
Executes the numerical integration to simulate the stochastic Wilson-Cowan dynamics.
- Parameters:
- pardict, optional
Dictionary of parameters to update before simulation
- x0array-like, optional
Initial state vector of length 2*nn (E and I for each region)
- verbosebool, default True
Whether to print simulation progress information
- Returns:
- tuple
Dictionary with keys ‘t’, ‘E’, ‘I’ containing time vector and excitatory/inhibitory time series arrays.
- integrate(P: ParWC, verbose=True)[source]¶
Run Wilson-Cowan model integration with Heun’s method.
Performs stochastic integration of the Wilson-Cowan equations using Heun’s method for improved accuracy. Returns time series data for excitatory and inhibitory populations.
- Parameters:
- PParWC
Parameter object containing all model parameters and settings
- verbosebool, default True
Whether to print integration progress
- Returns:
- dict
Dictionary containing: - ‘t’ : ndarray, time points - ‘E’ : ndarray or None, excitatory population activity - ‘I’ : ndarray or None, inhibitory population activity
- WC_sde¶
alias of
WC_sde_numba
vbi.models.numba.ww¶
- initialize_random_state(seed)[source]¶
Initialize the random number generator with a specific seed.
This function sets the random seed in the Numba context to ensure reproducible stochastic simulations.
- Parameters:
- seedint
Random seed value for reproducibility
- check_vec_size_1d(x, nn)[source]¶
Return a 1D vector of size nn, broadcasting scalar if needed.
This utility function ensures that parameter inputs are properly formatted as arrays of the correct size for multi-region simulations.
- Parameters:
- xscalar or array-like
Input value(s) to be broadcast or validated
- nnint
Required array size (number of brain regions)
- Returns:
- np.ndarray
Array of shape (nn,) with input values properly broadcast
- class ParBold(*args, **kwargs)[source]¶
Parameter class for BOLD signal generation in the Wong-Wang model.
This Numba jitclass holds parameters for the hemodynamic response model that converts neural activity to BOLD signal. Based on the Balloon-Windkessel model for simulating fMRI BOLD responses.
- Parameters:
- kappafloat, default 0.65
Signal decay parameter
- gammafloat, default 0.41
Feedback regulation parameter
- taufloat, default 0.98
Hemodynamic transit time (s)
- alphafloat, default 0.32
Grubb’s vessel stiffness exponent
- epsilonfloat, default 0.34
Efficacy of oxygen extraction
- Eofloat, default 0.4
Oxygen extraction fraction at rest
- TEfloat, default 0.04
Echo time (s)
- vofloat, default 0.08
Resting venous volume fraction
- r0float, default 25.0
Slope parameter for intravascular signal
- theta0float, default 40.3
Frequency offset at the outer surface of magnetized vessels (Hz)
- t_minfloat, default 0.0
Minimum integration time
- rtolfloat, default 1e-5
Relative tolerance for integration
- atolfloat, default 1e-8
Absolute tolerance for integration
- do_bold_step(r_in, s, f, ftilde, vtilde, qtilde, v, q, dtt, P)[source]¶
Perform one integration step of the BOLD hemodynamic response model.
Implements the Balloon-Windkessel model to convert neural activity into BOLD signal by simulating the hemodynamic cascade: neural activity → vasodilatory signal → blood flow → blood volume → deoxyhemoglobin → BOLD signal.
- Parameters:
- r_innp.ndarray
Neural activity input (excitatory synaptic gating variables)
- snp.ndarray, shape (2, nn)
Vasodilatory signal state variables
- fnp.ndarray, shape (2, nn)
Blood flow state variables
- ftildenp.ndarray, shape (2, nn)
Log-transformed blood flow variables
- vtildenp.ndarray, shape (2, nn)
Log-transformed blood volume variables
- qtildenp.ndarray, shape (2, nn)
Log-transformed deoxyhemoglobin content variables
- vnp.ndarray, shape (2, nn)
Blood volume state variables
- qnp.ndarray, shape (2, nn)
Deoxyhemoglobin content state variables
- dttfloat
Integration time step (seconds)
- PParBold
BOLD parameter object
Notes
This function modifies the state arrays in-place and uses a log-transform approach to ensure numerical stability of the hemodynamic state variables.
- class ParWW(*args, **kwargs)[source]¶
Parameter class for the Wong-Wang full neural mass model.
This Numba jitclass holds all parameters required for Wong-Wang simulation including biophysical parameters for excitatory and inhibitory populations, coupling strengths, noise levels, and simulation settings. Uses Numba for high-performance compilation.
- Parameters:
- a_excfloat, default 310.0
Excitatory population gain parameter (n/C)
- a_inhfloat, default 0.615
Inhibitory population gain parameter (nC⁻¹)
- b_excfloat, default 125.0
Excitatory population threshold parameter (Hz)
- b_inhfloat, default 177.0
Inhibitory population threshold parameter (Hz)
- d_excfloat, default 0.16
Excitatory population saturation parameter (s)
- d_inhfloat, default 0.087
Inhibitory population saturation parameter (s)
- tau_excfloat, default 100.0
Excitatory synaptic time constant (ms)
- tau_inhfloat, default 10.0
Inhibitory synaptic time constant (ms)
- gamma_excfloat, default 0.641/1000
Excitatory kinetic parameter (ms⁻¹)
- gamma_inhfloat, default 1.0/1000
Inhibitory kinetic parameter (ms⁻¹)
- W_excfloat, default 1.0
Excitatory population local weight
- W_inhfloat, default 0.7
Inhibitory population local weight
- ext_currentarray-like, default [0.382]
External current input per region (nA)
- J_NMDAfloat, default 0.15
NMDA synaptic coupling strength (nA)
- J_Ifloat, default 1.0
Inhibitory synaptic coupling strength (nA)
- w_plusfloat, default 1.4
Local excitatory recurrence strength
- lambda_inh_excfloat, default 0.0
Long-range feedforward inhibition switch
- firing_rate(current, a, b, d)[source]¶
Compute the firing rate using the Wong-Wang transfer function.
Implements the input-output relationship for neural populations: r(I) = (a*I - b) / (1 - exp(-d*(a*I - b)))
This function captures the nonlinear relationship between synaptic current and population firing rate, including saturation effects.
- Parameters:
- currentnp.ndarray
Synaptic current input to the population.
- afloat
Gain parameter controlling the slope of the transfer function.
- bfloat
Threshold parameter determining the firing threshold.
- dfloat
Saturation parameter controlling the saturation behavior.
- Returns:
- np.ndarray
Population firing rate (Hz).
- f_ww(S, t, P)[source]¶
Compute the right-hand side of the Wong-Wang full model equations.
This function implements the deterministic part of the Wong-Wang model, computing the time derivatives of synaptic gating variables for both excitatory and inhibitory populations across all brain regions.
The implemented equations are: dS_exc/dt = -S_exc/τ_exc + (1-S_exc)*γ_exc*r_exc(I_exc) dS_inh/dt = -S_inh/τ_inh + γ_inh*r_inh(I_inh)
- Parameters:
- Snp.ndarray
Current state vector of shape (2*nn,) containing stacked arrays: [S_exc_0, …, S_exc_n, S_inh_0, …, S_inh_n] where nn is the number of brain regions.
- tfloat
Current time (not used in autonomous system).
- PParWW
Parameter object containing all model parameters.
- Returns:
- np.ndarray
Derivative vector dS/dt of shape (2*nn,).
- heun_sde(S, t, P)[source]¶
Perform one heun integration step for the stochastic Wong-Wang model.
This function implements the Heun method (predictor-corrector) for numerical integration of stochastic differential equations, providing second-order accuracy for the deterministic part while properly handling additive Gaussian noise.
- Parameters:
- Snp.ndarray
Current state vector of shape (2*nn,) containing synaptic gating variables.
- tfloat
Current time (ms).
- PParWW
Parameter object containing model parameters including dt and sigma.
- Returns:
- np.ndarray
Updated state vector after one integration step.
- class WW_sde(par: dict | None = None, Bpar: dict | None = None)[source]¶
Numba implementation of the Wong-Wang full neural mass model with stochastic dynamics.
The Wong-Wang full model is a biophysically realistic neural mass model that explicitly captures the dynamics of both excitatory and inhibitory neural populations. This model is based on the original work of Wong and Wang [Wong2006] and has been extended for whole-brain network simulations [Deco2013], [Deco2014]. The model provides a detailed representation of recurrent network mechanisms underlying decision-making processes and has been widely used to study brain dynamics in health and disease.
The model describes the temporal evolution of synaptic gating variables for excitatory (S_exc) and inhibitory (S_inh) populations at each brain region. The dynamics are governed by the balance between synaptic decay, activity-dependent facilitation, and network coupling. The firing rates of each population are determined by input-output transfer functions that capture the relationship between synaptic currents and population firing rates.
The Wong-Wang full model equations are:
\[\begin{split}\frac{dS_{exc,i}}{dt} &= -\frac{S_{exc,i}}{\tau_{exc}} + (1 - S_{exc,i}) \gamma_{exc} r_{exc,i}(t) + \sigma \xi_i(t) \\ \frac{dS_{inh,i}}{dt} &= -\frac{S_{inh,i}}{\tau_{inh}} + \gamma_{inh} r_{inh,i}(t) + \sigma \xi_i(t)\end{split}\]where the firing rates are computed using:
\[\begin{split}r_{exc,i}(t) &= \frac{a_{exc} I_{exc,i} - b_{exc}}{1 - \exp(-d_{exc}(a_{exc} I_{exc,i} - b_{exc}))} \\ r_{inh,i}(t) &= \frac{a_{inh} I_{inh,i} - b_{inh}}{1 - \exp(-d_{inh}(a_{inh} I_{inh,i} - b_{inh}))}\end{split}\]The total synaptic currents for each population are:
\[\begin{split}I_{exc,i} &= W_{exc} I_{ext} + w_{plus} J_{NMDA} S_{exc,i} + G_{exc} J_{NMDA} \sum_{j=1}^{N} SC_{ij} S_{exc,j} - J_I S_{inh,i} \\ I_{inh,i} &= W_{inh} I_{ext} + J_{NMDA} S_{exc,i} - S_{inh,i} + G_{inh} J_{NMDA} \lambda_{inh,exc} \sum_{j=1}^{N} SC_{ij} S_{inh,j}\end{split}\]
Parameters¶ Name
Explanation
Default Value
a_excExcitatory population gain parameter (n/C)
310.0
a_inhInhibitory population gain parameter (nC⁻¹)
0.615
b_excExcitatory population threshold parameter (Hz)
125.0
b_inhInhibitory population threshold parameter (Hz)
177.0
d_excExcitatory population saturation parameter (s)
0.16
d_inhInhibitory population saturation parameter (s)
0.087
tau_excExcitatory synaptic time constant (ms)
100.0
tau_inhInhibitory synaptic time constant (ms)
10.0
gamma_excExcitatory kinetic parameter (ms⁻¹)
0.641/1000
gamma_inhInhibitory kinetic parameter (ms⁻¹)
1.0/1000
W_excExcitatory population local weight
1.0
W_inhInhibitory population local weight
0.7
ext_currentExternal current input (nA)
0.382
J_NMDANMDA synaptic coupling strength (nA)
0.15
J_IInhibitory synaptic coupling strength (nA)
1.0
w_plusLocal excitatory recurrence strength
1.4
lambda_inh_excLong-range feedforward inhibition switch
0.0
G_excGlobal excitatory coupling strength
0.0
G_inhGlobal inhibitory coupling strength
0.0
sigmaNoise amplitude
0.0
weightsStructural connectivity matrix (nn x nn)
zeros
dtIntegration time step (ms)
0.1
t_endSimulation end time (ms)
1000.0
t_cutInitial time to discard (ms)
0.0
- Usage example:
Notes
The Wong-Wang model is particularly useful for studying: - Decision-making processes and attractor dynamics - Excitatory-inhibitory balance in cortical circuits - Biophysically realistic neural mass dynamics - Working memory mechanisms - Brain state transitions and criticality
The model’s strength lies in its biophysical realism while maintaining computational efficiency through mean-field approximations.
References
[Wong2006]Wong, K. F., & Wang, X. J. (2006). A recurrent network mechanism of time integration in perceptual decisions. Journal of Neuroscience, 26(4), 1314-1328.
[Deco2013]Deco, G., Ponce-Alvarez, A., Mantini, D., Romani, G. L., Hagmann, P., & Corbetta, M. (2013). Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. Journal of Neuroscience, 33(27), 11239-11252.
[Deco2014]Deco, G., Ponce-Alvarez, A., Hagmann, P., Romani, G. L., Mantini, D., & Corbetta, M. (2014). How local excitation–inhibition ratio impacts the whole brain dynamics. Journal of Neuroscience, 34(23), 7886-7898.
- check_parameters(par: dict) None[source]¶
Validate that all provided parameters are recognized.
- Parameters:
- pardict
Dictionary of parameters to validate.
- Raises:
- ValueError
If any parameter name is not recognized.
- run(par: dict | None = None, x0=None, verbose=True)[source]¶
Run the Wong-Wang full model simulation.
Executes the Euler-Maruyama numerical integration scheme to simulate the stochastic Wong-Wang dynamics. Includes optional BOLD signal generation and data downsampling.
- Parameters:
- pardict, optional
Dictionary of parameters to update before simulation. Can include any parameter from the ParWW class.
- x0array-like, optional
Initial state vector of length 2*nn (S_exc and S_inh for each region). If None, uses the stored initial_state from the parameter object.
- verbosebool, default True
Whether to print simulation progress and parameter information.
- Returns:
- dict
Dictionary containing simulation results with keys:
- ‘S’ndarray, shape (T, nn)
Recorded excitatory synaptic gating variables if RECORD_S=True
- ‘t’ndarray, shape (T,)
Time points for neural data (ms)
- ‘bold_t’ndarray, shape (T_bold,)
Time points for BOLD signal (ms)
- ‘bold_d’ndarray, shape (T_bold, nn)
BOLD signal time series if RECORD_BOLD=True
Notes
The simulation uses the Heun-Euler method for stochastic integration, which provides second-order accuracy for the deterministic part while handling the stochastic noise appropriately.
- set_initial_state(nn, seed=-1)[source]¶
Generate random initial conditions for the Wong-Wang model.
Creates small positive random values for both excitatory and inhibitory synaptic gating variables, ensuring the system starts in a biologically plausible state.
- Parameters:
- nnint
Number of brain regions/nodes.
- seedint, optional
Random seed for reproducibility. If -1 or None, no seeding is applied.
- Returns:
- np.ndarray
Initial state vector of shape (2*nn,) containing small positive random values for [S_exc_0, …, S_exc_n, S_inh_0, …, S_inh_n].
PyTorch models¶
vbi.models.pytorch.ww_sde_kong¶
TVB-K models¶
vbi.models.tvbk.tvbk_wrapper¶
vbi.models.tvbk.utils¶
Core VBI modules¶
vbi.inference¶
vbi.utils¶
- timer(func)[source]¶
Decorator to measure elapsed time.
- Parameters:
- funcfunction
Function to be decorated.
- Returns:
- function
Wrapped function that measures execution time.
- display_time(time, message='')[source]¶
Display elapsed time in hours, minutes, seconds.
- Parameters:
- timefloat
Elapsed time in seconds.
- messagestr, optional
Optional message to display with the time. Default is empty string.
- class LoadSample(nn=84)[source]¶
Utility class for loading sample datasets and connectivity matrices.
This class provides convenient methods to load structural connectivity matrices, tract lengths, and BOLD signal data from the VBI dataset directory.
- Parameters:
- nnint, optional
Number of nodes/regions in the connectivity matrix. Default is 84. Supported values are typically 84 and 88.
- __init__(nn=84) None[source]¶
Initialize the LoadSample utility.
- Parameters:
- nnint, optional
Number of nodes/regions in the connectivity matrix. Default is 84.
- get_weights(normalize=True)[source]¶
Load structural connectivity weights matrix.
- Parameters:
- normalizebool, optional
Whether to normalize the weights by the maximum value. Default is True.
- Returns:
- np.ndarray
Structural connectivity matrix of shape (nn, nn) with diagonal set to 0. Values are non-negative after removing negative entries.
- get_limits(samples, limits=None)[source]¶
Calculate or validate parameter limits for samples.
This function computes the min/max limits for each parameter dimension across one or more sample arrays, or validates provided limits.
- Parameters:
- samplesnp.ndarray or list of np.ndarray
Sample array(s) of shape (n_samples, n_params) or list of such arrays. If PyTorch tensors, they will be converted to numpy arrays.
- limitslist or None, optional
Predefined limits as [[min1, max1], [min2, max2], …] for each parameter. If None or empty list, limits are computed from the data. If single limit pair provided, it will be broadcast to all parameters.
- Returns:
- torch.Tensor
Tensor of shape (n_params, 2) containing [min, max] for each parameter.
- get_limits_numpy(samples, limits=None)[source]¶
Calculate or validate parameter limits for samples (numpy-only version).
This function computes the min/max limits for each parameter dimension across one or more sample arrays, or validates provided limits.
- Parameters:
- samplesnp.ndarray or list of np.ndarray
Sample array(s) of shape (n_samples, n_params) or list of such arrays.
- limitslist or None, optional
Predefined limits as [[min1, max1], [min2, max2], …] for each parameter. If None or empty list, limits are computed from the data. If single limit pair provided, it will be broadcast to all parameters.
- Returns:
- np.ndarray
Array of shape (n_params, 2) containing [min, max] for each parameter.
- posterior_peaks(samples, return_dict=False, **kwargs)[source]¶
Find the peaks (modes) of a posterior distribution using kernel density estimation.
This function estimates the probability density of the posterior samples and identifies the locations of peak density for each parameter dimension.
- Parameters:
- samplesnp.ndarray or torch.Tensor
Posterior samples of shape (n_samples, n_params). If torch.Tensor, it will be converted to numpy array.
- return_dictbool, optional
If True, returns results as a dictionary with parameter labels as keys. If False, returns a simple list of peak values. Default is False.
- **kwargs
Additional keyword arguments passed to the plotting/analysis functions. These may include ‘labels’ for parameter names.
- Returns:
- list or dict
If return_dict=False: List of peak values for each parameter. If return_dict=True: Dictionary with parameter labels as keys and peak values as values.
- posterior_peaks_numpy(samples, return_dict=False, labels=None, bins=100, bw_method=None)[source]¶
Find the peaks (modes) of a posterior distribution using kernel density estimation (numpy-only version).
This function estimates the probability density of the posterior samples and identifies the locations of peak density for each parameter dimension. Does not require torch or sbi dependencies.
- Parameters:
- samplesnp.ndarray
Posterior samples of shape (n_samples, n_params).
- return_dictbool, optional
If True, returns results as a dictionary with parameter labels as keys. If False, returns a simple list of peak values. Default is False.
- labelslist or None, optional
Parameter names for dictionary keys. If None, uses integer indices.
- binsint, optional
Number of bins for density estimation grid. Default is 100.
- bw_methodstr, float or None, optional
Bandwidth method for KDE. Can be ‘scott’, ‘silverman’, or a scalar. If None, uses ‘scott’. Default is None.
- Returns:
- list or dict
If return_dict=False: List of peak values for each parameter. If return_dict=True: Dictionary with parameter labels as keys and peak values as values.
- Raises:
- ValueError
If samples contain insufficient data for KDE (need at least 2 samples).
- j2p(notebookPath, modulePath=None)[source]¶
convert a jupyter notebook to a python module
>>> j2p("sample.ipynb", "sample.py")
- posterior_shrinkage(prior_samples: torch.Tensor | np.ndarray, post_samples: torch.Tensor | np.ndarray) torch.Tensor[source]¶
Calculate the posterior shrinkage, quantifying how much the posterior distribution contracts from the initial prior distribution. References: https://arxiv.org/abs/1803.08393
- Parameters:
- prior_samplesarray_like or torch.Tensor [n_samples, n_params]
Samples from the prior distribution.
- post_samplesarray-like or torch.Tensor [n_samples, n_params]
Samples from the posterior distribution.
- Returns:
- shrinkagetorch.Tensor [n_params]
The posterior shrinkage.
- Raises:
- ImportError
If PyTorch is not installed.
- posterior_shrinkage_numpy(prior_samples: ndarray, post_samples: ndarray) ndarray[source]¶
Calculate the posterior shrinkage using numpy, quantifying how much the posterior distribution contracts from the initial prior distribution. References: https://arxiv.org/abs/1803.08393
- Parameters:
- prior_samplesnp.ndarray [n_samples, n_params]
Samples from the prior distribution.
- post_samplesnp.ndarray [n_samples, n_params]
Samples from the posterior distribution.
- Returns:
- shrinkagenp.ndarray [n_params]
The posterior shrinkage.
- Raises:
- ValueError
If input samples are empty or have incompatible shapes.
- posterior_zscore(true_theta: torch.Tensor | np.ndarray | float, post_samples: torch.Tensor | np.ndarray)[source]¶
Calculate the posterior z-score, quantifying how much the posterior distribution of a parameter encompasses its true value. References: https://arxiv.org/abs/1803.08393
- Parameters:
- true_thetafloat, array-like or torch.Tensor [n_params]
The true value of the parameters.
- post_samplesarray-like or torch.Tensor [n_samples, n_params]
Samples from the posterior distributions.
- Returns:
- ztorch.Tensor [n_params]
The z-score of the posterior distributions.
- Raises:
- ImportError
If PyTorch is not installed.
- posterior_zscore_numpy(true_theta: ndarray | float, post_samples: ndarray) ndarray[source]¶
Calculate the posterior z-score using numpy, quantifying how much the posterior distribution of a parameter encompasses its true value. References: https://arxiv.org/abs/1803.08393
- Parameters:
- true_thetafloat or np.ndarray [n_params]
The true value of the parameters.
- post_samplesnp.ndarray [n_samples, n_params]
Samples from the posterior distributions.
- Returns:
- znp.ndarray [n_params]
The z-score of the posterior distributions.
- Raises:
- ValueError
If input samples are empty.
- set_diag(A: ndarray, k: int = 0, value: float = 0.0)[source]¶
set k diagonals of the given matrix to given value.
- Parameters:
- A: np.ndarray
matrix
- k: int
number of diagonals
- value: float
value to be set
- Returns:
- A: np.ndarray
matrix with k diagonals set to value
- test_imports()[source]¶
Check required dependencies, including C++ modules, print versions, and warn if unavailable.
- get_cuda_info()[source]¶
Get CUDA version and device information using CuPy.
- Returns:
dict: Dictionary containing CUDA version and device information
- class BoxUniform(low, high, dtype=<class 'numpy.float64'>, seed=None)[source]¶
A multivariate uniform distribution over a hyperrectangle (box).
This implementation uses only numpy and scipy, providing a torch-free alternative for uniform distributions over rectangular domains.
- Parameters:
- lowarray_like
Lower bounds for each dimension. Shape (n_dims,) or scalar.
- higharray_like
Upper bounds for each dimension. Shape (n_dims,) or scalar.
- dtypenp.dtype, optional
Data type for internal arrays. Default is np.float64. Supported options are np.float32 and np.float64.
- seedint or None, optional
Random seed for reproducible sampling. If None, no seed is set.
- Attributes:
- lownp.ndarray
Lower bounds array.
- highnp.ndarray
Upper bounds array.
- ndimsint
Number of dimensions.
- volumefloat
Volume of the hyperrectangle.
- dtypenp.dtype
Data type used for arrays.
- __init__(low, high, dtype=<class 'numpy.float64'>, seed=None)[source]¶
Initialize BoxUniform distribution.
- Parameters:
- lowarray_like
Lower bounds for each dimension.
- higharray_like
Upper bounds for each dimension.
- dtypenp.dtype, optional
Data type for internal arrays. Default is np.float64. Supported options are np.float32 and np.float64.
- seedint or None, optional
Random seed for reproducible sampling. If None, no seed is set.
- set_seed(seed)[source]¶
Set the random seed for reproducible sampling.
- Parameters:
- seedint or None
Random seed. If None, the random state is not seeded.
- sample(sample_shape=(1,), seed=None)[source]¶
Sample from the uniform distribution.
- Parameters:
- sample_shapetuple or int, optional
Shape of samples to generate. Default is (1,).
- seedint or None, optional
Random seed for this sampling operation. If provided, it temporarily sets the seed for this operation only. If None, uses the current random state.
- Returns:
- np.ndarray
Samples of shape (*sample_shape, n_dims) with specified dtype.
- log_prob(value)[source]¶
Calculate log probability density.
- Parameters:
- valuearray_like
Values to evaluate. Shape (…, n_dims).
- Returns:
- np.ndarray
Log probability densities. Shape matches input except last dimension.
- prob(value)[source]¶
Calculate probability density.
- Parameters:
- valuearray_like
Values to evaluate. Shape (…, n_dims).
- Returns:
- np.ndarray
Probability densities. Shape matches input except last dimension.
- cdf(value)[source]¶
Calculate cumulative distribution function.
- Parameters:
- valuearray_like
Values to evaluate. Shape (…, n_dims).
- Returns:
- np.ndarray
CDF values. Shape matches input except last dimension.
- mean()[source]¶
Calculate the mean of the distribution.
- Returns:
- np.ndarray
Mean vector of shape (n_dims,).
- variance()[source]¶
Calculate the variance of the distribution.
- Returns:
- np.ndarray
Variance vector of shape (n_dims,).
vbi.optional_deps¶
Optional dependency handling for VBI.
This module provides utilities for gracefully handling optional dependencies and provides informative error messages when they’re missing.
- exception OptionalDependencyError[source]¶
Raised when an optional dependency is required but not available.
- optional_import(module_name: str, install_name: str | None = None) Any[source]¶
Import a module if available, otherwise return None.
- Parameters:
- module_namestr
Name of the module to import
- install_namestr, optional
Name to use in installation instructions (if different from module_name)
- Returns:
- module or None
The imported module if successful, None if not available
- require_optional(module_name: str, install_name: str | None = None, extra: str | None = None) Any[source]¶
Import a required optional dependency with helpful error message.
- Parameters:
- module_namestr
Name of the module to import
- install_namestr, optional
Name to use in installation instructions
- extrastr, optional
VBI extra that provides this dependency
- Returns:
- module
The imported module
- Raises:
- OptionalDependencyError
If the module cannot be imported
- requires_optional(*dependencies)[source]¶
Decorator to check for optional dependencies before function execution.
- Parameters:
- *dependenciestuples
Each tuple should be (module_name, install_name, extra)
Examples
>>> @requires_optional(('torch', 'torch', 'inference')) ... def inference_function(): ... import torch ... # function implementation ... pass
Feature extraction¶
vbi.feature_extraction.calc_features¶
- calc_features(ts: ndarray, fs: float, cfg: dict, preprocess=None, preprocess_args: dict = {}, verbose: bool = False)[source]¶
Extract features from time series data.
- Parameters:
- tsnp.ndarray
Time series data
- fsint, float
Sampling frequency
- cfgdict
Dictionary of features configurations
- preprocessfunction
Function for preprocessing the time series
- preprocess_argsdictionary
Arguments for preprocessing function
- Returns:
- featureslist of numpy arrays
- extract_features(ts: ~numpy.ndarray, fs: float, cfg: dict, output_type=<class 'vbi.feature_extraction.features_settings.Data_F'>, **kwargs)[source]¶
Extract features from time series data
- Parameters:
- tslist of np.ndarray [[n_regions x n_samples]]
Input from which the features are extracted
- fsint, float
Sampling frequency
- cfgdictionary
Dictionary of features to extract
- output_formatstring
Output format, either ‘list’ (list of numpy arrays) ‘dataframe’ (pandas dataframe) (default is ‘list’)
- **kwargs
- - n_workersint
Number of workers for parallelization, default is 1 Parallelization is done by ensembles (first dimension of ts)
- - dtypetype
Data type of the features extracted, default is np.float32
- - verboseboolean
If True, print the some information
- - preprocessfunction
Function for preprocessing the time series
- - preprocess_argsdictionary
Arguments for preprocessing function
- - n_trial: int
Number of trials
- Returns:
- Data Object with the following attributes:
- values: list of numpy arrays or pandas dataframe
extracted features
- labels: list of strings
List of labels of the features
- info: dictionary
Dictionary with the information of the features extracted
- extract_features_df(ts: ndarray, fs: float, cfg: dict, **kwargs)[source]¶
Extract features from time series data and return a pandas dataframe
- Parameters:
- tslist of np.ndarray [[n_regions x n_samples]]
Input from which the features are extracted
- fsint, float
Sampling frequency
- cfgdictionary
Dictionary of features to extract
- **kwargs
- - n_workersint
Number of workers for parallelization, default is 1 Parallelization is done by ensembles (first dimension of ts)
- - dtypetype
Data type of the features extracted, default is np.float32
- - verboseboolean
If True, print the some information
- - preprocessfunction
Function for preprocessing the time series
- - preprocess_argsdictionary
Arguments for preprocessing function
- Returns:
- Data Object with the following attributes:
- values: pandas dataframe
extracted features
- labels: list of strings
List of labels of the features
- info: dictionary
Dictionary with the information of the features extracted
- extract_features_list(ts, fs, cfg, **kwargs)[source]¶
extract features from time series data and return a list of features and labels
- Parameters:
- tslist of np.ndarray [[n_regions x n_samples]]
Input from which the features are extracted
- fsint, float
Sampling frequency
- cfgdictionary
Dictionary of features to extract
- **kwargs
- - n_workersint
Number of workers for parallelization, default is 1 Parallelization is done by ensembles (first dimension of ts)
- - dtypetype
Data type of the features extracted, default is np.float32
- - verboseboolean
If True, print the some information
- - preprocessfunction
Function for preprocessing the time series
- - preprocess_argsdictionary
Arguments for preprocessing function
- Returns:
- Data Object with the following attributes:
- values: list of numpy arrays
extracted features
- labels: list of strings
List of labels of the features
- info: dictionary
Dictionary with the information of the features extracted
vbi.feature_extraction.features¶
- abs_energy(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the absolute energy of the time serie.
>>> abs_energy([1, 2, 3, 4, 5]) (array([55]), ['abs_energy_0'])
- Parameters:
- tsnd-arrays [n_regions x n_samples]
Input from which the area under the curve is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: list of float
Absolute energy
- labels: list of str
Labels of the features
- average_power(ts: ndarray, fs: float = 1.0, indices: List[int] | None = None, verbose=False)[source]¶
Computes the average power of the time serie.
>>> average_power([1, 2, 3, 4, 5], 1) (array([13.75]), ['average_power_0'])
- Parameters:
- tsnd-arrays [n_regions x n_samples]
Input from which the area under the curve is computed
- fsfloat
Sampling frequency
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: list of float
Average power
- labels: list of str
Labels of the features
- auc(ts: ndarray, dx: float | None = None, x: ndarray | None = None, indices: List[int] | None = None, verbose=False)[source]¶
Computes the area under the curve of the signal computed with trapezoid rule.
>>> auc(np.array([[1, 2, 3], [4, 5, 6]]), None, np.array([0, 1, 2])) (array([ 4., 10.]), ['auc_0', 'auc_1'])
- Parameters:
- tsnd-arrays [n_regions x n_samples]
Input from which the area under the curve is computed
- dx: float
Spacing between values
- x: array_like, optional
x values of the time series
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- list of float
The area under the curve value
- labels: list of str
Labels of the features
- auc_lim(ts: ndarray, dx: float | None = None, x: ndarray | None = None, xlim: List[Tuple[float, float]] | None = None, indices: List[int] | None = None, verbose=False)[source]¶
Compute the area under the curve for a given time series within a given limit
>>> auc_lim(np.array([[1, 2, 3], [4, 5, 6]]), None, None, [(0, 1), (1, 2)]) ([1.5, 4.5, 2.5, 5.5], ['auc_lim_0', 'auc_lim_1', 'auc_lim_2', 'auc_lim_3'])
- Parameters:
- tsnd-arrays [n_regions x n_samples]
Input from which the area under the curve is computed
- dx: float
Spacing between values
- x: array_like
x values of the time series
- xlim: list of tuples
The limits of the time series
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- list of float
The area under the curve value
- labels: list of str
Labels of the features
- calc_var(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes variance of the time series.
>>> calc_var(np.array([[1, 2, 3], [4, 5, 6]])) (array([0.66666667, 0.66666667]), ['var_0', 'var_1'])
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which var is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
variance of the time series
- labels: array-like
labels of the features
- calc_std(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes standard deviation of the time serie.
>>> calc_std(np.array([[1, 2, 3], [4, 5, 6]])) (array([0.81649658, 0.81649658]), ['std_0', 'std_1'])
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which std is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
std of the time series
- labels: array-like
labels of the features
- calc_mean(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes median of the time serie.
>>> calc_mean(np.array([[1, 2, 3], [4, 5, 6]])) (array([2., 5.]), ['mean_0', 'mean_1'])
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which median is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
mean of the time series
- labels: array-like
labels of the features
- calc_centroid(ts: ndarray, fs: float, indices: List[int] | None = None, verbose=False)[source]¶
Computes the centroid along the time axis.
- Parameters:
- signalnd-array
Input from which centroid is computed
- fs: int
Signal sampling frequency
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- float
Temporal centroid
- calc_kurtosis(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the kurtosis of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which kurtosis is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
kurtosis of the time series
- labels: array-like
labels of the features
- calc_skewness(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the skewness of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which skewness is computed
- Returns:
- values: array-like
skewness of the time series
- labels: array-like
labels of the features
- calc_max(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the maximum of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which maximum is computed
- Returns:
- values: array-like
maximum of the time series
- labels: array-like
labels of the features
- calc_min(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the minimum of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which minimum is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
minimum of the time series
- labels: array-like
labels of the features
- calc_median(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the median of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which median is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
median of the time series
- labels: array-like
labels of the features
- mean_abs_dev(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the mean absolute deviation of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which mean absolute deviation is computed
- Returns:
- values: array-like
mean absolute deviation of the time series
- labels: array-like
labels of the features
- median_abs_dev(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the median absolute deviation of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which median absolute deviation is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
median absolute deviation of the time series
- labels: array-like
labels of the features
- rms(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the root mean square of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which root mean square is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
root mean square of the time series
- labels: array-like
labels of the features
- interq_range(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the interquartile range of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which interquartile range is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
interquartile range of the time series
- labels: array-like
labels of the features
- zero_crossing(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
Computes the number of zero crossings of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which number of zero crossings is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
number of zero crossings of the time series
- labels: array-like
labels of the features
- seizure_onset(ts: ndarray, threshold: float = 0.02, indices: List[int] | None = None, verbose=False)[source]¶
Computes the seizure onset of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which number of zero crossings is computed
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- values: array-like
index of the onset of the seizures in the time series, zero if no onset in each region
- labels: array-like
labels of the features
- kop(ts: ndarray, indices: List[int] | None = None, verbose=False, extract_phase=False)[source]¶
Calculate the Kuramoto order parameter (KOP)
The Kuramoto order parameter measures the synchronization level in a system of coupled oscillators. Values close to 1 indicate high synchronization, while values close to 0 indicate low synchronization.
- Parameters:
- tsnp.ndarray [n_regions x n_samples]
Input time series data
- indicesList[int], optional
Indices of the time series to compute the feature
- verbosebool, optional
Whether to print error messages
- extract_phasebool, optional
If True, extract phase information using Hilbert transform before computing KOP
- Returns:
- valueslist of float
Kuramoto order parameter values
- labelslist of str
Labels of the features
- calc_moments(ts: ndarray, indices: List[int] | None = None, orders: List[int] = [2, 3, 4, 5, 6], verbose=False)[source]¶
Computes the moments of the time series.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which moments are computed
- orders: list
List of orders of the moments
- Returns:
- values: array-like
moments of the time series
- labels: array-like
labels of the features
- calc_envelope(ts: ndarray, indices: List[int] | None = None, features: List[str] = ['mean', 'std', 'median', 'max', 'min'], verbose=False)[source]¶
Calculate statistics on the envelope of time series using Hilbert transform.
This function computes the analytic signal using Hilbert transform and extracts statistics from both the amplitude envelope and instantaneous phase.
- Parameters:
- tsnp.ndarray [n_regions x n_samples]
Input time series data
- indicesList[int], optional
Indices of the time series to compute the feature
- featuresList[str], optional
List of statistical features to compute on envelope Options: [“mean”, “std”, “median”, “max”, “min”]
- verbosebool, optional
Whether to print error messages
- Returns:
- valuesarray-like
Computed envelope statistics
- labelsarray-like
Labels of the features
- fc_sum(x: ndarray, positive=False, masks: Dict[str, ndarray] | None = None, verbose=False)[source]¶
Calculate the sum of functional connectivity (FC)
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which var is computed
- Returns:
- result: float
sum of functional connectivity
- fc_stat(ts: ndarray, k: int = 0, positive: bool = False, eigenvalues: bool = True, pca_num_components: int = 3, fc_function: str = 'corrcoef', masks: Dict[str, ndarray] | None = None, quantiles: List[float] = [0.05, 0.25, 0.5, 0.75, 0.95], features: List[str] = ['sum', 'max', 'min', 'mean', 'std', 'skew', 'kurtosis'], verbose=False)[source]¶
extract features from functional connectivity (FC)
- Parameters:
- ts: np.ndarray [n_regions, n_samples]
input array
- k: int
to remove up to kth diagonal of FC matrix
- pca_num_components: int
number of components for PCA
- positive: bool
if True, ignore negative values of fc elements
- masks: dict
dictionary of masks
- features: list of str
list of features to be extracted
- quantiles: list of float
list of quantiles, set 0 to ignore
- eigenvalues: bool
if True, extract features from eigenvalues
- fc_function: str
functional connectivity function: ‘corrcoef’ or ‘cov’
- Returns:
- stats: np.ndarray (1d)
feature values
- labels: list of str
feature labels
- fc_homotopic(ts: ndarray, average: bool = False, positive: bool = True, fc_function='corrcoef', verbose=False)[source]¶
Calculate the homotopic connectivity vector of a given brain activity
- Parameters:
- bold: array_like [nn, nt]
The brain activity to be analyzed.
- averag: bool
If True, the average homotopic connectivity is returned. Otherwise, the homotopic connectivity vector is returned.
- positive: bool
If True, only positive correlations are considered.
- Returns:
- valuesarray_like [n_nodes]
The homotopic correlation vector.
- labelslist of str
The labels of the homotopic correlation vector.
- Negative correlations may be artificially induced when using global signal regression
- in functional imaging pre-processing (Fox et al., 2009; Murphy et al., 2009; Murphy and Fox, 2017).
- Therefore, results on negative weights should be interpreted with caution and should be understood
- as complementary information underpinning the findings based on positive connections
- coactivation_degree(ts: ndarray, modality='noncor')[source]¶
Calculate coactivation degree (CAD). #! TODO need testing
Coactivation degree measures the temporal co-fluctuation of brain regions by computing the instantaneous product of regional activity with a global signal.
- Parameters:
- tsnp.ndarray [n_regions, n_samples]
Input time series array
- modalitystr, optional
Modality for global signal computation - “noncor”: Exclude current region from global signal (default) - “cor”: Include all regions in global signal
- Returns:
- valueslist
Coactivation degree values for each region-timepoint pair
- labelslist
Labels of the features (empty list as this returns raw values)
Notes
This function is currently under development and testing.
- coactivation_phase(ts)[source]¶
Calculate the coactivation phase (CAP). # ! TODO need testing
Coactivation phase measures the phase relationship between regional signals and the global signal using Hilbert transform to extract instantaneous phases.
- Parameters:
- tsnp.ndarray [n_regions, n_samples]
Input time series array
- Returns:
- CAPlist
Mean phase differences between regional and global signals
Notes
This function is currently under development and testing. The function computes instantaneous phases using Hilbert transform and calculates the mean phase difference for each region.
- burstiness(ts: ndarray, indices: List[int] | None = None, verbose=False)[source]¶
calculate the burstiness statistic - Goh and Barabasi, ‘Burstiness and memory in complex systems’ Europhys. Lett. 81, 48002 (2008). [from hctsa-py]
- Parameters:
- x: np.ndarray [n_regions, n_samples]
input array
- indices: list of int
Indices of the time series to compute the feature
- Returns:
- B: list of floats
burstiness statistic
- fcd_stat(ts, TR=1, win_len=30, masks=None, positive=False, eigenvalues=True, pca_num_components=3, quantiles=[0.05, 0.25, 0.5, 0.75, 0.95], features=['sum', 'max', 'min', 'mean', 'std', 'skew', 'kurtosis'], k=None, verbose=False)[source]¶
- calc_mi(ts: ndarray, k: int = 4, time_diff: int = 1, num_threads: int = 1, source_indices: List[int] | None = None, target_indices: List[int] | None = None, mode: str = 'pairwise', verbose=False, **kwargs)[source]¶
calculate the mutual information between time series based on the Kraskov method #!TODO bug in multiprocessing
- Parameters:
- ts: np.ndarray [n_regions, n_samples]
input array
- k: int
kth nearest neighbor
- time_diff: int
time difference between time series
- num_threads: int
number of threads
- source_indices: list or np.ndarray
indices of source time series, if None, all time series are used
- target_indices: list or np.ndarray
indices of target time series, if None, all time series are used
- mode: str
“pairwise” or “all”, if “pairwise”, source_indices and target_indices must have the same length
- Returns:
- MI: list of floats
mutual information
- labels: list of str
labels of the features
- calc_te(ts: ndarray, k: int = 4, delay: int = 1, num_threads: int = 1, source_indices: List[int] | None = None, target_indices: List[int] | None = None, mode: str = 'pairwise', verbose=False, **kwargs)[source]¶
calculate the transfer entropy between time series based on the Kraskov method.
- Parameters:
- ts: np.ndarray [n_regions, n_samples]
input array
- num_threads: int
number of threads
- source_indices: list or np.ndarray
indices of source time series, if None, all time series are used
- target_indices: list or np.ndarray
indices of target time series, if None, all time series are used
- mode: str
“pairwise” or “all”, if “pairwise”, source_indices and target_indices must have the same length
- Returns:
- TE: list of floats
transfer entropy
- calc_entropy(ts: ndarray, average: bool = False, verbose=False)[source]¶
Calculate entropy of time series using Kozachenko-Leonenko estimator.
This function computes the differential entropy of the time series data using the Kozachenko-Leonenko k-nearest neighbor entropy estimator.
- Parameters:
- tsnp.ndarray [n_regions x n_samples]
Input time series data
- averagebool, optional
If True, compute average entropy across all regions If False, compute entropy for each region separately
- verbosebool, optional
Whether to print error messages
- Returns:
- valueslist of float or float
Entropy values in bits
- labelslist of str or str
Labels of the features
- calc_entropy_bin(ts: ndarray, prob: str = 'standard', average: bool = False, verbose=False)[source]¶
Computes the entropy of the signal using the Shannon Entropy.
Description in Article: Regularities Unseen, Randomness Observed: Levels of Entropy Convergence Authors: Crutchfield J. Feldman David
- Parameters:
- signalnd-array
Input from which entropy is computed
- probstring
Probability function (kde or gaussian functions are available)
- Returns:
- values: float or array-like
The normalized entropy value
- labels: string or array-like
The label of the feature
- spectrum_stats(ts: ndarray, fs: float, method: str = 'fft', nperseg: int | None = None, verbose=False, indices: List[int] | None = None, average=False, features: List[str] = ['spectral_distance', 'fundamental_frequency', 'max_frequency', 'max_psd', 'median_frequency', 'spectral_centroid', 'spectral_kurtosis', 'spectral_variation'])[source]¶
Compute various statistics of the power spectrum of time series.
This function calculates multiple spectral features including spectral distance, fundamental frequency, maximum frequency, maximum PSD, median frequency, spectral centroid, spectral kurtosis, and spectral variation.
- Parameters:
- tsnp.ndarray [n_regions x n_samples]
Input time series from which power spectrum statistics are computed
- fsfloat
Sampling frequency in Hz
- methodstr, optional
Method to compute the power spectrum. Options: ‘welch’, ‘fft’ (default: ‘fft’)
- npersegint, optional
Length of each segment for Welch method. If None, uses half the time series length
- verbosebool, optional
Whether to print error messages
- indicesList[int], optional
Indices of the regions to be used. If None, all regions are used
- averagebool, optional
If True, average PSD across regions before computing features
- featuresList[str], optional
List of spectral features to compute
- Returns:
- valuesarray-like
Computed power spectrum statistics
- labelsarray-like
Labels of the features
- spectrum_auc(ts, fs, method='fft', bands=None, nperseg=None, average=False, indices=None, verbose=False)[source]¶
calculate the area under the curve of the power spectrum of the time series over given frequency bands.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input time series
- fsfloat
Sampling frequency
- methodstr
Method to compute the power spectrum. Can be ‘welch’ or ‘fft’
- bandslist of tuples
Frequency bands
- nperseg: int
Length of each segment. default is half of the time series
- avg: bool
averaging psd over all regions
- indices: list of int
indices of the regions to be used
- Returns:
- values: array-like
area under the curve of the power spectrum of the time series
- labels: array-like
labels of the features
- spectrum_moments(ts, fs, method='fft', nperseg=None, avg=False, moments=[2, 3, 4, 5, 6], normalize=False, indices=None, average=False, verbose=False)[source]¶
Computes the moments of power spectrum
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which power spectrum statistics are computed
- fsfloat
Sampling frequency
- methodstr
Method to compute the power spectrum. Can be ‘welch’ or ‘fft’
- nperseg: int
…
- avg: bool
averaging over all regions
- nm: list of int
moments orders
- Returns:
- values: array-like
power spectrum statistics of the time series
- labels: array-like
labels of the features
- psd_raw(ts, fs, bands=[(0, 4), (4, 8), (8, 12), (12, 30), (30, 70)], df=None, method='fft', nperseg=None, average=False, normalize=False, normalize_to: float | None = None, indices=None, verbose=False)[source]¶
Calculate frequency spectrum and return with specified frequency resolution.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input time series
- fsfloat
Sampling frequency
- bandslist of tuples
Frequency bands
- dffloat
Frequency resolution, default is fs / n_samples
- methodstr
Method to compute the power spectrum. Can be ‘welch’ or ‘fft’
- nperseg: int
Length of each segment. default is half of the time series
- avg: bool
averaging psd over all regions
- normalize: bool
normalize the psd by the maximum value
- normalize_to: float
normalize the psd to the given frequency value
- indices: list of int
indices of the regions to be used
- Returns:
- psd: array-like
power spectrum density
- wavelet_abs_mean_1d(ts, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]), verbose=False)[source]¶
Computes CWT absolute mean value of each wavelet scale.
- Parameters:
- tsnd-array
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- tuple
CWT absolute mean value
- wavelet_abs_mean(ts, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]), verbose=False)[source]¶
Computes CWT absolute mean value of each wavelet scale.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- values: array-like
CWT absolute mean value of the time series
- labels: array-like
labels of the features
- wavelet_std(ts, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]), verbose=False)[source]¶
Computes CWT std value of each wavelet scale.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- values: array-like
CWT std value of the time series
- labels: array-like
labels of the features
- wavelet_energy_1d(ts, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]), verbose=False)[source]¶
Computes CWT energy of each wavelet scale.
Implementation details: https://stackoverflow.com/questions/37659422/energy-for-1-d-wavelet-in-python
Feature computational cost: 2
- Parameters:
- signalnd-array
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- tuple
CWT energy
- wavelet_energy(ts, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]), verbose=False)[source]¶
Computes CWT energy of each wavelet scale.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- values: array-like
CWT energy of the time series
- labels: array-like
labels of the features
- hmm_stat(ts, node_indices=None, n_states=4, subname='', n_iter=100, seed=None, observations='gaussian', method='em', tcut=5, bins=10, verbose=False)[source]¶
Calculate Hidden Markov Model (HMM) statistics including state durations and transition matrix.
This function fits an HMM to the time series data and extracts features related to state durations and transition probabilities.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input time series from which HMM features are computed
- node_indiceslist, optional
List of node indices to be used for HMM fitting If None, all nodes are used
- n_statesint, optional
Number of hidden states (default: 4)
- subnamestr, optional
Substring to add to feature labels
- n_iterint, optional
Number of EM iterations for fitting (default: 100)
- seedint, optional
Random seed for reproducibility
- observationsstr, optional
Observation distribution type (default: “gaussian”)
- methodstr, optional
Method to fit the HMM (default: “em”)
- tcutint, optional
Maximum duration of a state for histogram (default: 5)
- binsint, optional
Number of bins for state duration histogram (default: 10)
- verbosebool, optional
Whether to print verbose output
- Returns:
- stat_vecarray-like
Concatenated HMM features (state durations + transition matrix)
- labelsarray-like
Labels of the features
- catch22(ts, indices: List[int] | None = None, catch24=False, verbose=False, features=['DN_HistogramMode_5', 'DN_HistogramMode_10', 'CO_f1ecac', 'CO_FirstMin_ac', 'CO_HistogramAMI_even_2_5', 'CO_trev_1_num', 'MD_hrv_classic_pnn40', 'SB_BinaryStats_mean_longstretch1', 'SB_TransitionMatrix_3ac_sumdiagcov', 'PD_PeriodicityWang_th0_01', 'CO_Embed2_Dist_tau_d_expfit_meandiff', 'IN_AutoMutualInfoStats_40_gaussian_fmmi', 'FC_LocalSimple_mean1_tauresrat', 'DN_OutlierInclude_p_001_mdrmd', 'DN_OutlierInclude_n_001_mdrmd', 'SP_Summaries_welch_rect_area_5_1', 'SB_BinaryStats_diff_longstretch0', 'SB_MotifThree_quantile_hh', 'SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1', 'SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1', 'SP_Summaries_welch_rect_centroid', 'FC_LocalSimple_mean3_stderr'])[source]¶
Calculate the Catch22 features.
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which Catch22 features are computed
- node_indiceslist
List of node indices to be used for Catch22
- catch24bool
If True, calculate mean and std of the features
- Returns:
- valuesarray-like
feature values
- labelsarray-like
labels of the features
vbi.feature_extraction.features_settings¶
- load_json(path)[source]¶
Load json file
- Parameters:
- pathstring
Path to json file
- Returns:
- json_datadictionary
Dictionary with the json data
- get_features_by_domain(domain=None, json_path=None)[source]¶
Create a dictionary of features in given domain(s).
- Parameters:
- domainlist of strings or None, optional
List of domains of features to extract. If None, all domains are returned. Valid domains are: ‘hmm’, ‘spectral’, ‘connectivity’, ‘temporal’, ‘statistical’, ‘information’, ‘catch22’.
- json_pathstring or None, optional
Path to json file containing feature definitions. If None, uses the default features.json file in the package.
- Returns:
- dict
Dictionary of features filtered by the specified domain(s). Keys are domain names, values are dictionaries of features in that domain.
- get_features_by_given_names(cfg, names=None)[source]¶
Filter features by given names from cfg (a dictionary of features).
- Parameters:
- cfgdict
Dictionary of features organized by domain. Each domain contains features as key-value pairs.
- nameslist of strings, tuple of strings, string, or None, optional
Names of features to extract. Can be a single name (string) or multiple names (list/tuple). If None, returns the original cfg. Names are case-insensitive.
- Returns:
- dict
Dictionary of features filtered by the specified names. Structure matches input cfg but only contains features with matching names.
Notes
If a feature name is not found in the available features, a warning message is printed but processing continues.
- get_features_by_tag(tag=None, json_path=None)[source]¶
Create a dictionary of features in given tag.
- Parameters:
- tagstring or None, optional
Tag of features to extract. Valid tags include “fmri”, “audio”, “eeg”, “ecg”. If None, returns all features.
- json_pathstring or None, optional
Path to json file containing feature definitions. If None, uses the default features.json file in the package.
- Returns:
- dict
Dictionary of features filtered by the specified tag. Keys are domain names, values are dictionaries of features in that domain that match the tag. Empty domains are removed from the result.
- Raises:
- SystemExit
If tag is not one of the valid options: “audio”, “eeg”, “ecg”, or None.
- add_feature(cfg, domain, name, function: str | None = None, features_path: str | ModuleType | None = None, parameters={}, tag=None, description='')[source]¶
Add a feature to the cfg dictionary
- Parameters:
- cfgdictionary
Dictionary of features
- domainstring
Domain of the feature
- namestring
Name of the feature
- functionfunction
Function to compute the feature
- parametersdictionary
Parameters of the feature
- tagstring
Tag of the feature
- descriptionstring
Description of the feature
- Returns:
- cfgdictionary
Updated dictionary of features
- add_features_from_json(json_path, features_path, fea_dict={})[source]¶
Add features from json file to cfg dictionary.
- Parameters:
- json_pathstring
Path to json file containing feature definitions to load.
- features_pathstring or module
Path to the module containing the feature functions, or the module object itself.
- fea_dictdict, optional
Dictionary of features to add to. If empty, a new dictionary is created. Default is {}.
- Returns:
- dict
Dictionary containing all features from the json file added to the input fea_dict.
Notes
TODO: Check if features already exist in fea_dict to avoid conflicts. TODO: Check for conflicts in parameters and function definitions.
- class Data_F(values=None, labels=None, info=None)[source]¶
Data container class for feature extraction results.
A simple container to store feature values along with their labels and additional information.
- Parameters:
- valuesarray-like or None, optional
Feature values, typically a numpy array or list. Default is None.
- labelsarray-like or None, optional
Labels corresponding to the feature values. Default is None.
- infodict or None, optional
Additional information about the features (e.g., parameters used, feature names, etc.). Default is None.
- Attributes:
- valuesarray-like or None
The feature values.
- labelsarray-like or None
The feature labels.
- infodict or None
Additional feature information.
- update_cfg(cfg: dict, name: str, parameters: dict)[source]¶
Set parameters of a feature in the configuration dictionary.
- Parameters:
- cfgdict
Dictionary of features organized by domain. Each domain contains features as key-value pairs.
- namestr
Name of the feature to update.
- parametersdict
Parameters as key-value pairs to set for the feature.
- Returns:
- dict
Updated dictionary of features with the specified feature’s parameters modified.
Notes
This function searches through all domains to find the feature by name and updates its parameters. If the feature is not found, the cfg is returned unchanged.
- select_features_by_domain(module_name, domain)[source]¶
Select functions from a module that belong to a specific domain.
Note: This function is not currently used in the codebase.
- Parameters:
- module_namestr
Name of the module to inspect for functions.
- domainstr
Domain name to filter functions by. Functions must have a ‘domain’ attribute containing this value.
- Returns:
- list
List of function objects that have the specified domain in their ‘domain’ attribute.
- select_functions_by_tag(module_name, tag)[source]¶
Select functions from a module that have a specific tag.
Note: This function is not currently used in the codebase.
- Parameters:
- module_namestr
Name of the module to inspect for functions.
- tagstr
Tag name to filter functions by. Functions must have a ‘tag’ attribute containing this value.
- Returns:
- list
List of function objects that have the specified tag in their ‘tag’ attribute.
- select_functions_by_domain_and_tag(module_name, domain=None, tag=None)[source]¶
Select functions from a module that match both domain and tag criteria.
Note: This function is not currently used in the codebase.
- Parameters:
- module_namestr
Name of the module to inspect for functions.
- domainstr or None, optional
Domain name to filter functions by. Functions must have a ‘domain’ attribute containing this value.
- tagstr or None, optional
Tag name to filter functions by. Functions must have a ‘tag’ attribute containing this value.
- Returns:
- list
List of function objects that have both the specified domain in their ‘domain’ attribute and the specified tag in their ‘tag’ attribute.
vbi.feature_extraction.features_utils¶
- slice_features(x: ndarray | None, feature_names: list, info: dict)[source]¶
Slice features using given feature list
- Parameters:
- x: array-like
- features: list of strings
list of features
- info: dict
features’s colum indices in x
- Returns:
- x_sliced: array-like
sliced features
- preprocess(ts, fs=None, preprocess_dict={}, **kwargs)[source]¶
Preprocess time series data
- Parameters:
- tsnd-array [n_regions, n_timepoints]
Input from which the features are extracted
- fsint
Sampling frequency, set to 1 if not used
- preprocess_dictdictionary
Dictionary of preprocessing options
- **kwargsdict
Additional arguments
- band_pass_filter(ts, low_cut=0.02, high_cut=0.1, TR=2.0, order=2)[source]¶
apply band pass filter to given time series
- Parameters:
- tsnumpy.ndarray [n_regions, n_timepoints]
Input signal
- low_cutfloat, optional
Low cut frequency. The default is 0.02.
- high_cutfloat, optional
High cut frequency. The default is 0.1.
- TRfloat, optional
Sampling interval. The default is 2.0 second.
- Returns:
- ts_filtnumpy.ndarray
filtered signal
- remove_strong_artefacts(ts, threshold=3.0)[source]¶
Remove strong artifacts from time series by clipping values beyond threshold.
This function identifies outlier values in the time series that exceed a certain number of standard deviations and clips them to the threshold value to reduce the impact of artifacts on subsequent analysis.
- Parameters:
- tsarray-like [n_regions, n_timepoints] or list
Input time series data
- thresholdfloat, optional
Number of standard deviations beyond which values are considered artifacts Default is 3.0
- Returns:
- tsnp.ndarray [n_regions, n_timepoints]
Time series with artifacts removed (clipped to threshold)
Examples
>>> import numpy as np >>> ts = np.array([[1, 2, 100, 4, 5], [2, 3, 4, -50, 6]]) >>> clean_ts = remove_strong_artefacts(ts, threshold=2.0)
- get_fc(ts, masks=None, positive=False, fc_fucntion='corrcoef')[source]¶
calculate the functional connectivity matrix
- Parameters:
- tsnumpy.ndarray [n_regions, n_timepoints]
Input signal
- Returns:
- FCnumpy.ndarray
functional connectivity matrix
- get_fcd(ts, TR=1, win_len=30, positive=False, masks=None)[source]¶
Compute dynamic functional connectivity.
- Parameters:
- ts: numpy.ndarray [n_regions, n_timepoints]
Input signal
- win_len: int
sliding window length in samples, default is 30
- TR: int
repetition time. It refers to the amount of time that passes between consecutive acquired brain volumes during functional magnetic resonance imaging (fMRI) scans.
- positive: bool
if True, only positive values of FC are considered. default is False
- masks: dict
dictionary of masks to compute FCD on. default is None, which means that FCD is computed on the full matrix. see also
hbt.utility.make_maskandhbt.utility.get_masks.- Returns:
- FCD: ndarray
matrix of functional connectivity dynamics
- get_fcd2(ts, wwidth=30, maxNwindows=200, olap=0.94, indices=[], verbose=False)[source]¶
Calculate Functional Connectivity Dynamics (FCD) from time series using sliding windows.
This function computes the dynamic functional connectivity by calculating correlation matrices in sliding time windows and then computing the correlation between these windowed connectivity patterns over time.
- Parameters:
- tsnp.ndarray [n_nodes, n_samples]
Input time series data with nodes as rows and time samples as columns
- wwidthint, optional
Window width in time samples (default: 30)
- maxNwindowsint, optional
Maximum number of windows to compute (default: 200)
- olapfloat, optional
Overlap between consecutive windows as fraction (0-1, default: 0.94)
- indiceslist, optional
List of node indices to include in analysis (default: empty list uses all)
- verbosebool, optional
Whether to print verbose output (default: False)
- Returns:
- FCDnp.ndarray [n_windows, n_windows]
Functional connectivity dynamics matrix representing correlations between windowed connectivity patterns
Notes
The FCD matrix captures how functional connectivity patterns change over time by correlating the upper triangular elements of windowed FC matrices.
- filterbank(signal, fs, pre_emphasis=0.97, nfft=512, nfilt=40)[source]¶
Computes the MEL-spaced filterbank.
It provides the information about the power in each frequency band.
Implementation details and description on: https://www.kaggle.com/ilyamich/mfcc-implementation-and-tutorial https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html#fnref:1
- Parameters:
- signalnd-array
Input from which filterbank is computed
- fsint
Sampling frequency
- pre_emphasisfloat
Pre-emphasis coefficient for pre-emphasis filter application
- nfftint
Number of points of fft
- nfiltint
Number of filters
- Returns:
- nd-array
MEL-spaced filterbank
- autocorr_norm(signal)[source]¶
Compute normalized autocorrelation function of a signal.
This function calculates the autocorrelation of a signal normalized by the variance and length, providing a measure of how similar the signal is to shifted versions of itself.
- Parameters:
- signalnp.ndarray
Input signal from which autocorrelation is computed. Should be a 1D array.
- Returns:
- acfnp.ndarray
Normalized autocorrelation function of the same length as input signal. Values range from 0 to 1, where 1 indicates perfect correlation at lag 0.
Notes
Implementation details and description in: https://ccrma.stanford.edu/~orchi/Documents/speaker_recognition_report.pdf
The autocorrelation is normalized by variance and signal length to provide a standardized measure independent of signal amplitude and duration.
Examples
>>> import numpy as np >>> signal = np.sin(np.linspace(0, 4*np.pi, 100)) >>> acf = autocorr_norm(signal) >>> print(acf[0]) # Should be close to 1.0
- create_symmetric_matrix(acf, order=11)[source]¶
Computes a symmetric matrix.
Implementation details and description in: https://ccrma.stanford.edu/~orchi/Documents/speaker_recognition_report.pdf
- Parameters:
- acfnd-array
Input from which a symmetric matrix is computed
- orderint
Order
- Returns:
- nd-array
Symmetric Matrix
- lpc(signal, n_coeff=12)[source]¶
Computes the linear prediction coefficients.
Implementation details and description in: https://ccrma.stanford.edu/~orchi/Documents/speaker_recognition_report.pdf
- Parameters:
- signalnd-array
Input from linear prediction coefficients are computed
- n_coeffint
Number of coefficients
- Returns:
- nd-array
Linear prediction coefficients
- create_xx(features)[source]¶
Computes the range of features amplitude for the probability density function calculus.
- Parameters:
- featuresnd-array
Input features
- Returns:
- nd-array
range of features amplitude
- kde(features)[source]¶
Compute probability density function using Gaussian Kernel Density Estimation.
This function estimates the probability density function of the input data using a Gaussian KDE with Silverman’s bandwidth selection method.
- Parameters:
- featuresnp.ndarray
Input data from which probability density function is computed. Should be a 1D array of numerical values.
- Returns:
- pdfnp.ndarray
Normalized probability density values corresponding to the input range. Sum of all values equals 1.
Notes
Uses Silverman’s rule-of-thumb for bandwidth selection
Adds small noise if all values are identical to avoid singularity
Evaluates PDF over linearly spaced points covering the data range
Examples
>>> import numpy as np >>> data = np.random.normal(0, 1, 100) >>> pdf = kde(data) >>> print(np.sum(pdf)) # Should be approximately 1.0
- gaussian(features)[source]¶
Compute probability density function using a fitted Gaussian distribution.
This function fits a Gaussian (normal) distribution to the input data and evaluates the probability density function over the data range.
- Parameters:
- featuresnp.ndarray
Input data from which probability density function is computed. Should be a 1D array of numerical values.
- Returns:
- pdfnp.ndarray
Normalized probability density values from the fitted Gaussian distribution. Sum of all values approximates 1.
Notes
Fits a normal distribution using sample mean and standard deviation
Evaluates PDF over linearly spaced points covering the data range
More parametric than KDE but assumes Gaussian underlying distribution
Examples
>>> import numpy as np >>> data = np.random.normal(5, 2, 100) >>> pdf = gaussian(data) >>> print(pdf.shape) # Same length as input data
- calc_ecdf(signal)[source]¶
- Computes the ECDF of the signal.
ECDF is the empirical cumulative distribution function.
- Parameters:
- signalnd-array
Input from which ECDF is computed
- Returns
- ——-
- nd-array
Sorted signal and computed ECDF.
- matrix_stat(A: ndarray, k: int = 1, eigenvalues: bool = True, pca_num_components: int = 3, quantiles: List[float] = [0.05, 0.25, 0.5, 0.75, 0.95], features: List[str] = ['sum', 'max', 'min', 'mean', 'std', 'skew', 'kurtosis'])[source]¶
Calculate comprehensive statistics from a matrix (typically connectivity matrices).
This function extracts various statistical features from a matrix including basic statistics, eigenvalue properties, PCA components, and quantiles. Commonly used for analyzing functional connectivity matrices.
- Parameters:
- Anp.ndarray [n x n]
Input matrix, typically a square connectivity or correlation matrix
- kint, optional
Upper triangular matrix offset. Only elements above the k-th diagonal are considered (default: 1, excludes main diagonal)
- eigenvaluesbool, optional
Whether to compute eigenvalue-based features (default: True)
- pca_num_componentsint, optional
Number of PCA components to extract. Set to 0 to skip PCA (default: 3)
- quantilesList[float], optional
List of quantiles to compute (default: [0.05, 0.25, 0.5, 0.75, 0.95]) Set to [] or None to skip quantile computation
- featuresList[str], optional
List of statistical features to compute from matrix values Options: [“sum”, “max”, “min”, “mean”, “std”, “skew”, “kurtosis”]
- Returns:
- valuesnp.ndarray
Concatenated array of all computed feature values
- labelslist of str
Corresponding feature labels describing each value
Notes
This function is particularly useful for connectivity analysis where you need to extract summary statistics from FC/SC matrices while avoiding redundant information from symmetric matrices.
Examples
>>> import numpy as np >>> # Create a sample correlation matrix >>> A = np.random.rand(10, 10) >>> A = (A + A.T) / 2 # Make symmetric >>> values, labels = matrix_stat(A, k=1) >>> print(f"Computed {len(values)} features")
- init_jvm()[source]¶
Initialize Java Virtual Machine for information theory calculations.
- Raises:
- ImportError
If JPype is not available
- compute_time(ts, fs)[source]¶
Creates the signal correspondent time array.
- Parameters:
- signal: nd-array
Input from which the time is computed.
- fs: int
Sampling Frequency
- Returns:
- timefloat list
Signal time
- calc_fft(ts, fs)[source]¶
This functions computes the fft of a signal.
- Parameters:
- signalnd-array [n_regions, n_timepoints]
The input signal from which fft is computed
- fsfloat
Sampling frequency
- Returns:
- f: nd-array
Frequency values (xx axis)
- fmag: nd-array [n_regions, n_freqs]
Amplitude of the frequency values (yy axis)
- fundamental_frequency(f, fmag)[source]¶
Computes fundamental frequency of the signal.
The fundamental frequency integer multiple best explain the content of the signal spectrum.
Feature computational cost: 1
- Parameters:
- tsnd-array [n_regions x n_samples]
Input from which fundamental frequency is computed
- fsfloat
Sampling frequency
- Returns:
- f0: array of floats
Predominant frequency of the signals
- spectral_distance(freq, fmag)[source]¶
Computes the signal spectral distance.
Distance of the signal’s cumulative sum of the FFT elements to the respective linear regression.
- Parameters:
- fmag: nd-array [n_regions x n_freqs]
power spectrum of the signal
- Returns:
- values: array-like
spectral distances
- labels: array-like
labels of the features
- max_frequency(f, psd)[source]¶
Computes the maximum frequency of the signals.
- Parameters:
- f: nd-array
frequency values
- psd: nd-array [n_regions x n_freqs]
power spectral density of the signal
- Returns:
- values: array-like
maximum frequencies
- max_psd(f, psd)[source]¶
Computes the maximum power spectral density of the signals.
- Parameters:
- f: nd-array
frequency values
- psd: nd-array [n_regions x n_freqs]
power spectral density of the signal
- Returns:
- values: array-like
maximum power spectral densities
- spectral_centroid(f, fmag)[source]¶
Calculate the spectral centroid of the signals. The Timbre Toolbox: Extracting audio descriptors from musicalsignals Authors Peeters G., Giordano B., Misdariis P., McAdams S.
- Parameters:
- f: nd-array
frequency values
- fmag: nd-array [n_regions x n_freqs]
power spectrum of the signal
- Returns:
- values: array-like
spectral centroids
- labels: array-like
labels of the features
- spectral_kurtosis(f, fmag)[source]¶
Measure the flatness of the power spectrum of the signals. The Timbre Toolbox: Extracting audio descriptors from musicalsignals Authors Peeters G., Giordano B., Misdariis P., McAdams S.
- Parameters:
- f: nd-array
frequency values
- fmag: nd-array [n_regions x n_freqs]
power spectrum of the signal
- Returns:
- values: array-like
spectral kurtosis
- labels: array-like
labels of the features
- spectral_spread(f, fmag)[source]¶
Measures the spread of the spectrum around its mean value.
Description and formula in Article: The Timbre Toolbox: Extracting audio descriptors from musicalsignals Authors Peeters G., Giordano B., Misdariis P., McAdams S.
Feature computational cost: 2
- Parameters:
- signalnd-array
Signal from which spectral spread is computed.
- fsfloat
Sampling frequency
- Returns:
- float
Spectral Spread
- spectral_variation(freq, fmag)[source]¶
Computes the amount of variation of the spectrum along time. Spectral variation is computed from the normalized cross-correlation between two consecutive amplitude spectra.
Description and formula in Article: The Timbre Toolbox: Extracting audio descriptors from musicalsignals Authors Peeters G., Giordano B., Misdariis P., McAdams S.
- wavelet(signal, function=None, widths=array([1, 2, 3, 4, 5, 6, 7, 8, 9]))[source]¶
Computes CWT (continuous wavelet transform) of the signal.
- Parameters:
- signalnd-array
Input from which CWT is computed
- functionwavelet function
Default: scipy.signal.ricker
- widthsnd-array
Widths to use for transformation Default: np.arange(1,10)
- Returns:
- nd-array
The result of the CWT along the time axis matrix with size (len(widths),len(signal))
- km_order(ts, indices=None, avg=True)[source]¶
Calculate the (local) Kuramoto order parameter (KOP) of the given time series
- Parameters:
- ts: np.ndarray (2d) [n_regions, n_timepoints]
input array
- indices: list
list of indices of the regions of interest
- avg: bool
if True, average the KOP across time
- Returns:
- values: np.ndarray (1d) or float
feature values
- normalize_signal(ts, method='zscore')[source]¶
Normalize the input time series
- Parameters:
- ts: np.ndarray (2d) [n_regions, n_timepoints]
input array
- method: str
normalization method
- index: int
index of the times point to normalize with respect to x = x / x[:, index]
- Returns:
- ts: np.ndarray (2d) [n_regions, n_timepoints]
normalized array
- state_duration(hmm_z: ndarray, n_states: int, avg: bool = True, tcut: int = 5, bins: int = 10)[source]¶
Measure the duration of each state
- Parameters:
- hmm_znd-array [n_samples]
The most likely states for each time point
- n_statesint
The number of states
- avgbool
If True, the average duration of each state is returned. Otherwise, the duration of each state is returned.
- t_cutint
maximum duration of a state, default is 5
- binsint
number of bins for the histogram, default is 10
- Returns:
- stat_vecarray-like
The duration of each state
vbi.feature_extraction.utility¶
- prepare_input(ts, dtype=<class 'numpy.float32'>)[source]¶
prepare input format
- Parameters:
- tsarray-like or list
Input from which the features are extracted
- Returns
- ——-
- ts: nd-array
formatted input
- make_mask(n, indices)[source]¶
make a mask matrix with given indices
- Parameters:
- nint
size of the mask matrix
- indiceslist
indices of the mask matrix
- Returns:
- masknumpy.ndarray
mask matrix
- get_intrah_mask(n_nodes)[source]¶
Get a mask for intrahemispheric connections.
- Parameters:
- n_nodes: int
number of total nodes that constitute the data.
- Returns:
- mask_intrah: 2d array
mask for intrahemispheric connections.
- get_interh_mask(n_nodes)[source]¶
Get a mask for interhemispheric connections.
- Parameters:
- n_nodes: int
number of total nodes that constitute the data.
- Returns:
- mask_interh: 2d array
mask for interhemispheric connections.
- get_masks(n_nodes, networks)[source]¶
Get a dictionary of masks based on the requested networks.
- Parameters:
- n_nodes: int
number of total nodes that constitute the data.
- networks: list of str
list of networks to be included in the dictionary. ‘full’: full-network connections ‘intrah’: intrahemispheric connections ‘interh’: interhemispheric connections to get a custom mask with specific indices refere to
hbt.utility.make_mask(n, indices).- Returns:
- masks: dict
dictionary of masks based on the requested networks.
- is_sequence(arg)[source]¶
Check if the input is a sequence (list, tuple, np.ndarray, etc.)
- Parameters:
- argany
input to be checked.
- Returns:
- bool
True if the input is a sequence, False otherwise.
- set_k_diagonals(A, k=0, value=0)[source]¶
set k diagonals of the given matrix to given value.
- Parameters:
- Anumpy.ndarray
input matrix.
- kint
number of diagonals to be set. The default is 0. Notice that the main diagonal is 0.
- valueint, optional
value to be set. The default is 0.
- if_symmetric(A, tol=1e-08)[source]¶
Check if the input matrix is symmetric.
- Parameters:
- Anumpy.ndarray
input matrix.
- tolfloat, optional
tolerance for checking symmetry. The default is 1e-8.
- Returns:
- bool
True if the input matrix is symmetric, False otherwise.
- scipy_iir_filter_data(x, sfreq, l_freq, h_freq, l_trans_bandwidth=None, h_trans_bandwidth=None, **kwargs)[source]¶
Custom, scipy based filtering function with basic butterworth filter. #comes from neurolib
- Parameters:
- xnp.ndarray
data to be filtered, time is the last axis
- sfreqfloat
sampling frequency of the data in Hz
- l_freqfloat|None
frequency below which to filter the data in Hz
- h_freqfloat|None
frequency above which to filter the data in Hz
- l_trans_bandwidthkeeping for compatibility with mne
- h_trans_bandwidthkeeping for compatibility with mne
- **kwargspossible keywords to
scipy.signal.butter:- Returns:
- np.ndarray
filtered data
- filter(ts: ndarray, fs: float, low_freq: float, high_freq: float, l_trans_bandwidth: str = 'auto', h_trans_bandwidth: str = 'auto', **kwargs)[source]¶
- Filter data. Can be:
low-pass (low_freq is None, high_freq is not None),
high-pass (high_freq is None, low_freq is not None),
band-pass (l_freq < h_freq),
band-stop (l_freq > h_freq) filter type
- Parameters:
- ts: np.ndarray
Time series data
- low_freqfloat|None
frequency below which to filter the data.
- high_freqfloat|None
frequency above which to filter the data.
- l_trans_bandwidthfloat|str
transition band width for low frequency
- h_trans_bandwidthfloat|str
transition band width for high frequency
- inplacebool
whether to do the operation in place or return
- kwargspossible keywords to mne.filter.create_filter:
filter_length=”auto”, method=”fir”, iir_params=None phase=”zero”, fir_window=”hamming”, fir_design=”firwin”
- Returns:
- np.ndarray
filtered data
- posterior_shrinkage(prior_samples: None | ndarray, post_samples: None | ndarray) None[source]¶
Calculate the posterior shrinkage, quantifying how much the posterior distribution contracts from the initial prior distribution. References: https://arxiv.org/abs/1803.08393
- Parameters:
- prior_samplesarray_like or torch.Tensor [n_samples, n_params]
Samples from the prior distribution.
- post_samplesarray-like or torch.Tensor [n_samples, n_params]
Samples from the posterior distribution.
- Returns:
- shrinkagetorch.Tensor [n_params]
The posterior shrinkage.
- posterior_zscore(true_theta: None | array | float, post_samples: None | array)[source]¶
Calculate the posterior z-score, quantifying how much the posterior distribution of a parameter encompasses its true value. References: https://arxiv.org/abs/1803.08393
- Parameters:
- true_thetafloat, array-like or torch.Tensor [n_params]
The true value of the parameters.
- post_samplesarray-like or torch.Tensor [n_samples, n_params]
Samples from the posterior distributions.
- Returns:
- zTensor [n_params]
The z-score of the posterior distributions.