import warnings
import numpy as np
from numba import njit
from numba.experimental import jitclass
from numba.core.errors import NumbaPerformanceWarning
from numba import float64, int64
warnings.simplefilter("ignore", category=NumbaPerformanceWarning)
[docs]
@njit(nogil=True)
def run(P, times):
"""
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
----------
P : ParGHB
Parameter object containing all model parameters.
times : np.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).
"""
G = P.G
dt = P.dt
eta = P.eta
SC = P.weights
sigma = P.sigma
omega = P.omega
tcut = P.tcut
decimate = P.decimate
init_state = P.init_state
epsilon = 0.5
itaus = 1.25
itauf = 2.5
itauo = 1.02040816327
ialpha = 5.0
Eo = 0.4
V0 = 4.0
k1 = 2.77264
k2 = 0.572
k3 = -0.43
nt = times.shape[0]
nn = SC.shape[0]
# state variables
n_buffer = np.int64(np.floor(nt / decimate) + 1)
x = np.zeros(nn)
bold = np.zeros((nn, n_buffer))
y = np.zeros(nn)
z = np.array([0.0] * nn + [1.0] * 3 * nn)
# act = np.zeros((nn, nt))
# initial conditions (similar value for all regions)
x_init, y_init = init_state[:nn], init_state[nn:]
x[:] = x_init
y[:] = y_init
for i in range(nn):
bold[i, 0] = V0 * (
k1
- k1 * z[3 * nn + i]
+ k2
- k2 * (z[3 * nn + i] / z[2 * nn + i])
+ k3
- k3 * z[2 * nn + i]
)
ii = 0 # counter for decimation
for it in range(nt - 1):
for i in range(nn):
gx, gy = 0.0, 0.0
for j in range(nn):
gx = gx + SC[i, j] * (x[j] - x[i])
gy = gy + SC[i, j] * (y[j] - y[i])
dx = (
(x[i] * (eta[i] - (x[i] * x[i]) - (y[i] * y[i])))
- (omega[i] * y[i])
+ (G * gx)
)
dy = (
(y[i] * (eta[i] - (x[i] * x[i]) - (y[i] * y[i])))
+ (omega[i] * x[i])
+ (G * gy)
)
dz0 = epsilon * x[i] - itaus * z[i] - itauf * (z[nn + i] - 1)
dz1 = z[i]
dz2 = itauo * (z[nn + i] - z[2 * nn + i] ** ialpha)
dz3 = itauo * (
z[nn + i] * (1 - (1 - Eo) ** (1 / z[nn + i])) / Eo
- (z[2 * nn + i] ** ialpha) * z[3 * nn + i] / z[2 * nn + i]
)
x[i] = x[i] + dt * dx + np.sqrt(dt) * sigma * np.random.randn()
y[i] = y[i] + dt * dy + np.sqrt(dt) * sigma * np.random.randn()
z[i] = z[i] + dt * dz0
z[nn + i] = z[nn + i] + dt * dz1
z[2 * nn + i] = z[2 * nn + i] + dt * dz2
z[3 * nn + i] = z[3 * nn + i] + dt * dz3
if (it%decimate == 0):
bold[i, ii + 1] = V0 * (
k1
- k1 * z[3 * nn + i]
+ k2
- k2 * (z[3 * nn + i] / z[2 * nn + i])
+ k3
- k3 * z[2 * nn + i]
)
if (it%decimate == 0):
ii += 1
bold = bold[:, times[::decimate]>tcut]
t_bold = times[times[::decimate]>tcut]
return t_bold, bold
[docs]
class GHB_sde(object):
"""
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
.. list-table:: Parameters
:widths: 25 50 25
:header-rows: 1
* - Name
- Explanation
- Default Value
* - `G`
- Global coupling strength scaling network connections.
- 1.0
* - `dt`
- Integration time step in seconds.
- 0.001
* - `sigma`
- Noise amplitude (standard deviation of additive Gaussian noise).
- 0.1
* - `tend`
- End time of simulation in seconds.
- 10.0
* - `tcut`
- Time from which to start collecting output (burn-in period).
- 0.0
* - `eta`
- Bifurcation parameter array of length `nn`. Controls oscillation amplitude. If array-like, it should be of length `nn` (number of nodes).
- np.array([])
* - `omega`
- Intrinsic frequency array of length `nn`. Sets the natural oscillation frequency for each region. If array-like, it should be of length `nn`.
- np.array([])
* - `weights`
- Structural connectivity matrix of shape (`nn`, `nn`). Must be provided.
- np.array([[], []])
* - `decimate`
- Decimation factor for output time series (every `decimate`-th point is saved).
- 1
* - `init_state`
- Initial state vector of shape (2*nn,) containing [x₀, y₀]. If None, random initial conditions are generated.
- np.array([])
* - `seed`
- Random seed for reproducible simulations. If -1, no seeding is applied.
- -1
Usage example:
>>> import numpy as np
>>> from vbi.models.numba.ghb import GHB_sde
>>> W = np.eye(2) # 2-node demo connectivity
>>> eta = np.array([0.1, 0.1]) # Bifurcation parameters
>>> omega = np.array([40.0, 40.0]) # Frequencies in Hz
>>> ghb = GHB_sde({
... "weights": W,
... "eta": eta,
... "omega": omega,
... "dt": 0.01,
... "tend": 10.0,
... "tcut": 2.0
... })
>>> result = ghb.run()
>>> t, bold = result["t"], result["bold"] # bold has shape (nn, n_timepoints)
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
"""
[docs]
def __init__(self, par: dict = {}) -> None:
"""
Initialize the Generic Hopf Bifurcation model.
Parameters
----------
par : dict, optional
Dictionary containing model parameters. See class documentation for
available parameters. The 'weights', 'eta', and 'omega' parameters
are required for proper functionality.
"""
self.valid_par = [par_spec[i][0] for i in range(len(par_spec))]
self.check_parameters(par)
self.P = self.get_par_obj(par)
[docs]
def get_par_obj(self, par: dict):
"""
Create parameter object from dictionary.
Parameters
----------
par : dict
Parameter dictionary.
Returns
-------
ParGHB
Numba jitclass parameter object.
"""
if "init_state" in par.keys():
par["init_state"] = np.array(par["init_state"])
if "weights" in par.keys():
par["weights"] = np.array(par["weights"])
return ParGHB(**par)
def __str__(self) -> str:
"""
Return a string representation of the model parameters.
Returns
-------
str
Formatted string showing all model parameters and their values.
"""
print("Generic Hopf Bifurcation Model (Numba)")
print("---------------------------------------")
for key in self.valid_par:
print(f"{key}: {getattr(self.P, key)}")
return ""
[docs]
def check_parameters(self, par: dict) -> None:
"""
Validate that all provided parameters are recognized.
Parameters
----------
par : dict
Dictionary of parameters to validate.
Raises
------
ValueError
If any parameter name is not recognized.
"""
for key in par.keys():
if key not in self.valid_par:
raise ValueError(f"Invalid parameter: {key}")
[docs]
def set_initial_state(self, seed=None):
"""
Generate random initial conditions for the model.
Parameters
----------
seed : int, optional
Random seed for reproducible initial conditions.
Returns
-------
np.ndarray
Initial state vector of shape (2*nn,) with random values in [0, 1].
"""
if seed is not None:
np.random.seed(seed)
assert self.P.weights is not None
return np.random.uniform(0, 1, 2 * self.P.weights.shape[0])
[docs]
def run(self, par={}, tspan=None, x0=None, verbose=True):
"""
Run the Generic Hopf Bifurcation simulation.
Parameters
----------
par : dict, optional
Dictionary of parameters to update for this simulation run.
Any parameter from the class documentation can be updated.
tspan : tuple or None, optional
Time span as (start_time, end_time). If None, uses (0, tend).
x0 : np.ndarray, optional
Initial state vector of shape (2*nn,). If None, random initial
conditions are generated using set_initial_state().
verbose : bool, 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
"""
if x0 is None:
self.seed = self.P.seed if self.P.seed > 0 else None
self.P.init_state = self.set_initial_state(seed=self.seed)
else:
self.P.init_state = x0
if tspan is None:
times = np.arange(0, self.P.tend, self.P.dt)
else:
times = np.arange(tspan[0], tspan[1], self.P.dt)
if par:
self.check_parameters(par)
for key in par.keys():
setattr(self.P, key, par[key])
self.check_input()
t, b = run(self.P, times)
return {'t': t, 'bold': b}
par_spec = [
("G", float64),
("dt", float64),
("seed", int64),
("tend", float64),
("tcut", float64),
("sigma", float64),
("eta", float64[:]),
("decimate", int64),
("omega", float64[:]),
("weights", float64[:, :]),
("init_state", float64[:]),
]
[docs]
@jitclass(par_spec)
class ParGHB:
"""
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.
"""
def __init__(
self,
G=1.0,
dt=0.001,
sigma=0.1,
tend=10.0,
tcut=0.0,
eta=np.array([]),
init_state=np.array([]),
omega=np.array([]),
weights=np.array([[], []]),
decimate=1,
):
"""
Initialize the parameter container.
Parameters
----------
G : float
Global coupling strength.
dt : float
Integration time step.
sigma : float
Noise amplitude.
tend : float
End time of simulation.
tcut : float
Burn-in time.
eta : np.ndarray
Bifurcation parameters for each node.
init_state : np.ndarray
Initial state vector.
omega : np.ndarray
Intrinsic frequencies for each node.
weights : np.ndarray
Connectivity matrix.
decimate : int
Decimation factor.
"""
self.G = G
self.dt = dt
self.seed = -1
self.eta = eta
self.tend = tend
self.tcut = tcut
self.sigma = sigma
self.omega = omega
self.weights = weights
self.decimate = decimate
self.init_state = init_state
# Alias for consistency with naming convention
GHB_sde_numba = GHB_sde