PyRates Parameter Analysis

This guide demonstrates parameter sweeps and grid searches using TVBO’s SimulationExperiment.explorations with PyRates as the backend. Grid searches are declared declaratively — no manual loops, no YAML export plumbing. TVBO translates each Exploration into a PyRates-native grid_search() call automatically.

Workflow

  1. Define the model with Dynamics
  2. Declare the exploration with SimulationExperiment.explorations
  3. Run exp.run('pyrates') — TVBO invokes PyRates grid_search() natively
  4. Access results via result.explorations[name] (ExplorationResult)

Setup

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

Example 1: Van der Pol Parameter Sweep

Explore how the damping parameter μ affects oscillation behaviour.

Define Model and Exploration

from tvbo import Dynamics, SimulationExperiment
from tvbo.datamodel.schema import Exploration, Parameter, Range
from IPython.display import Markdown, display

# Van der Pol oscillator
vdp = Dynamics("Dynamics")
vdp.name = "VdP"
vdp.add_parameter("mu", value=1.0, description="Damping parameter")
vdp.add_state_variable("x", equation="z", initial_value=0.1)
vdp.add_state_variable("z", equation="mu*(1 - x**2)*z - x", initial_value=0.0)

display(Markdown(vdp.generate_report(format="markdown")))

VdP

State Equations

\[ \dot{x} = z \] \[ \dot{z} = - x + \mu*z*\left(1 - x^{2}\right) \]

Parameters

Parameter Value Unit Description
\(\mu\) 1.0 N/A Damping parameter

Declare Exploration and Run

exp = SimulationExperiment(dynamics=vdp)
exp.integration.duration = 50.0
exp.integration.step_size = 1e-2

# Declare 1-D sweep over mu — TVBO translates this to PyRates grid_search()
exp.explorations['mu_sweep'] = Exploration(
    name='mu_sweep',
    label='μ sweep',
    description='Effect of damping coefficient μ on VdP oscillations',
    parameters=[Parameter(name='mu', domain=Range(lo=0.1, hi=10.0, n=6))],
    mode='product',
)

print("Running μ sweep via PyRates grid_search …")
result = exp.run('pyrates')
sweep = result.explorations['mu_sweep']

mu_values = sweep.axes[0].explored_values
n_time = sweep.results.shape[1]
time = np.arange(n_time) * sweep.dt
x_idx = sweep.output_names.index('x')
z_idx = sweep.output_names.index('z')

print(f"Sweep complete — {len(mu_values)} conditions × {n_time} time steps")
print(f"ExplorationResult: {sweep}")
Running μ sweep via PyRates grid_search …
Compilation Progress
--------------------
    (1) Translating the circuit template into a networkx graph representation...
        ...finished.
    (2) Preprocessing edge transmission operations...
        ...finished.
    (3) Parsing the model equations into a compute graph...
        ...finished.
    Model compilation was finished.
Simulation Progress
-------------------
     (1) Generating the network run function...
     (2) Processing output variables...
        ...finished.
     (3) Running the simulation...
        ...finished after 0.02262170799986052s.
Compilation Progress
--------------------
    (1) Translating the circuit template into a networkx graph representation...
        ...finished.
    (2) Preprocessing edge transmission operations...
        ...finished.
    (3) Parsing the model equations into a compute graph...
        ...finished.
    Model compilation was finished.
Simulation Progress
-------------------
     (1) Generating the network run function...
     (2) Processing output variables...
        ...finished.
     (3) Running the simulation...
        ...finished after 0.05941066699961084s.
Sweep complete — 6 conditions × 5000 time steps
ExplorationResult: ExplorationResult(name='mu_sweep', grid=6, timeseries=(6, 5000, 2))

Visualize Time Series

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()

for i, mu in enumerate(mu_values):
    ax = axes[i]
    ax.plot(time, sweep.results[i, :, x_idx], label='x')
    ax.plot(time, sweep.results[i, :, z_idx], alpha=0.7, label='z')
    ax.set_title(f'μ = {mu:.2g}')
    ax.set_xlabel('Time')
    ax.set_ylabel('State')
    ax.legend(loc='upper right')

plt.suptitle('Van der Pol: Effect of Damping (μ)')
plt.tight_layout()
plt.show()

Phase Portraits

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
axes = axes.flatten()
colors = cm.viridis(np.linspace(0, 1, len(mu_values)))

for i, (mu, color) in enumerate(zip(mu_values, colors)):
    ax = axes[i]
    x = sweep.results[i, :, x_idx]
    z = sweep.results[i, :, z_idx]
    ax.plot(x, z, color=color, linewidth=0.8)
    ax.scatter(x[0], z[0], color='green', s=50, zorder=5, label='Start')
    ax.scatter(x[-1], z[-1], color='red', s=50, zorder=5, label='End')
    ax.set_title(f'μ = {mu:.2g}')
    ax.set_xlabel('x')
    ax.set_ylabel('z')
    ax.set_box_aspect(1)

plt.suptitle('Van der Pol Phase Portraits')
plt.tight_layout()
plt.show()

Example 3: Bifurcation-like Sweep

Scan the background drive η of a QIF mean-field model to find the transition from silence to oscillation.

Define Model and Sweep

qif = Dynamics("Dynamics")
qif.name = "QIF"
qif.add_parameter("tau",   value=1.0)
qif.add_parameter("eta",   value=-5.0,  description="Background drive")
qif.add_parameter("J",     value=15.0,  description="Coupling strength")
qif.add_parameter("Delta", value=1.0,   description="Heterogeneity")
qif.add_state_variable("r", equation="Delta/(tau*pi) + 2*r*v",         initial_value=0.1)
qif.add_state_variable("v", equation="v**2 + eta + J*r*tau - (pi*r*tau)**2", initial_value=-2.0)

display(Markdown(qif.generate_report(format="markdown")))

QIF

State Equations

\[ \dot{r} = 2*r*v + \frac{\Delta}{\pi*\tau} \] \[ \dot{v} = \eta + v^{2} + J*r*\tau - \pi^{2}*r^{2}*\tau^{2} \]

Parameters

Parameter Value Unit Description
\(\tau\) 1.0 N/A None
\(\eta\) -5.0 N/A Background drive
\(J\) 15.0 N/A Coupling strength
\(\Delta\) 1.0 N/A Heterogeneity 
exp3 = SimulationExperiment(dynamics=qif)
exp3.integration.duration = 200.0
exp3.integration.step_size = 5e-3

exp3.explorations['eta_sweep'] = Exploration(
    name='eta_sweep',
    label='η bifurcation sweep',
    parameters=[Parameter(name='eta', domain=Range(lo=-10.0, hi=5.0, n=15))],
    mode='product',
)

print("Running η bifurcation sweep …")
result3  = exp3.run('pyrates')
sweep3   = result3.explorations['eta_sweep']

eta_values = sweep3.axes[0].explored_values
n_time3    = sweep3.results.shape[1]
r_idx      = sweep3.output_names.index('r')

print(f"Sweep complete — {len(eta_values)} conditions × {n_time3} time steps")
Running η bifurcation sweep …
Compilation Progress
--------------------
    (1) Translating the circuit template into a networkx graph representation...
        ...finished.
    (2) Preprocessing edge transmission operations...
        ...finished.
    (3) Parsing the model equations into a compute graph...
        ...finished.
    Model compilation was finished.
Simulation Progress
-------------------
     (1) Generating the network run function...
     (2) Processing output variables...
        ...finished.
     (3) Running the simulation...
        ...finished after 0.385886207999647s.
Compilation Progress
--------------------
    (1) Translating the circuit template into a networkx graph representation...
        ...finished.
    (2) Preprocessing edge transmission operations...
        ...finished.
    (3) Parsing the model equations into a compute graph...
        ...finished.
    Model compilation was finished.
Simulation Progress
-------------------
     (1) Generating the network run function...
     (2) Processing output variables...
        ...finished.
     (3) Running the simulation...
        ...finished after 4.052057500000956s.
Sweep complete — 15 conditions × 40000 time steps

Bifurcation Diagram

steady_start = int(0.8 * n_time3)

r_min  = [float(np.min( sweep3.results[i, steady_start:, r_idx])) for i in range(len(eta_values))]
r_max  = [float(np.max( sweep3.results[i, steady_start:, r_idx])) for i in range(len(eta_values))]
r_mean = [float(np.mean(sweep3.results[i, steady_start:, r_idx])) for i in range(len(eta_values))]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
ax = axes[0]
ax.fill_between(eta_values, r_min, r_max, alpha=0.3, label='Oscillation range')
ax.plot(eta_values, r_mean, 'k-', linewidth=2, label='Mean firing rate')
ax.set_xlabel('η (background drive)')
ax.set_ylabel('Firing rate r')
ax.set_title('QIF: Bifurcation-like Diagram')
ax.legend()
ax.grid(True, alpha=0.3)

time3 = np.arange(n_time3) * sweep3.dt
ax = axes[1]
sample_indices = [0, 5, 9, 14]
colors_ = plt.cm.viridis(np.linspace(0, 1, len(sample_indices)))
for idx, color in zip(sample_indices, colors_):
    start = int(0.6 * n_time3)
    ax.plot(time3[start:] - time3[start], sweep3.results[idx, start:, r_idx],
            color=color, label=f'η = {eta_values[idx]:.1f}')
ax.set_xlabel('Time')
ax.set_ylabel('Firing rate r')
ax.set_title('Sample Time Series')
ax.legend()

plt.tight_layout()
plt.show()

Summary

TVBO’s explorations declaration maps directly to PyRates grid_search():

TVBO concept PyRates concept
Exploration.parameters + Range(lo, hi, n) param_grid = {name: linspace(lo, hi, n)}
Exploration.mode = 'product' permute_grid=True (Cartesian product)
Exploration.mode = 'zip' permute_grid=False (pairwise)
result.explorations[name] returned (results_df, results_map) DataFrames 
from tvbo import Dynamics, SimulationExperiment
from tvbo.datamodel.schema import Exploration, Parameter, Range

exp = SimulationExperiment(dynamics=model)
exp.integration.duration = T
exp.integration.step_size = dt

exp.explorations['my_sweep'] = Exploration(
    name='my_sweep',
    parameters=[Parameter(name='mu', domain=Range(lo=0.1, hi=10.0, n=6))],
    mode='product',  # or 'zip'
)

result = exp.run('pyrates')
sweep  = result.explorations['my_sweep']

# Access:  sweep.axes[i].explored_values, sweep.results[cond, time, var], sweep.output_names

Next Steps