Note
Go to the end to download the full example code.
Optimizing Descending Drive Parameters for Target Firing Rates#
This example demonstrates parameter optimization for descending drive networks using Optuna. The goal is to find network parameters (number of DD neurons, connection probability, drive frequency) that produce motor neuron firing patterns matching target physiological characteristics.
Note
Multi-objective optimization is essential for tuning complex neural network models because:
Multiple parameters interact non-linearly (DD neurons, connectivity, synaptic weights)
Target outputs have multiple constraints (mean firing rate, variability, recruitment)
Manual parameter tuning is time-consuming and may miss optimal combinations
Systematic search ensures reproducible, well-documented parameter choices
Important
Descending Drive (DD) Networks represent cortical input to spinal motor neurons. Key parameters:
dd_neurons: Population size (affects input diversity and convergence)conn_probability: Synaptic connectivity (affects drive strength and correlation)dd_drive__Hz: Input frequency (controls baseline excitation level)gamma_shape: Variability of Poisson processes (affects temporal patterns)
Optimization Objective: Match target motor neuron firing rate statistics while maintaining biologically plausible network parameters.
Import Libraries#
For optimization, we need:
Optuna: Bayesian optimization framework with TPE sampler
NEURON: Biophysical simulation engine
MyoGen: Motor neuron pools and network connectivity
import json
import os
os.environ["MPLBACKEND"] = "Agg"
if "DISPLAY" in os.environ:
del os.environ["DISPLAY"]
import warnings
from pathlib import Path
import joblib
import matplotlib.pyplot as plt
import numpy as np
import optuna
import quantities as pq
from neo import Segment, SpikeTrain
from neuron import h
from myogen import RANDOM_GENERATOR
from myogen.simulator import RecruitmentThresholds
from myogen.simulator.neuron import Network
from myogen.simulator.neuron.populations import AlphaMN__Pool, DescendingDrive__Pool
from myogen.utils.helper import calculate_firing_rate_statistics, get_gamma_shape_for_mvc
from myogen.utils.nmodl import load_nmodl_mechanisms
warnings.filterwarnings("ignore")
plt.style.use("fivethirtyeight")
Define Optimization Parameters#
Set target firing rate characteristics and optimization constraints. These values are based on physiological measurements from motor control studies.
# Target motor neuron firing statistics (from experimental data)
TARGET_FR_MEAN__HZ = 16.82 # Mean firing rate during isometric contraction
TARGET_FR_STD__HZ = 2.5 # Standard deviation across motor unit pool
# Network parameter search space
N_DD_NEURONS_MIN = 100 # Minimum descending drive population size
N_DD_NEURONS_MAX = 1000 # Maximum descending drive population size
# Fixed simulation parameters
N_MOTOR_UNITS = 100 # Motor neuron pool size
SIMULATION_TIME_MS = 3000.0 # Simulation duration (ms)
TIMESTEP_MS = 0.1 # Integration timestep (ms)
GAMMA_SHAPE = 3.0 # Fixed gamma shape for DD variability
SYNAPTIC_WEIGHT = 0.05 # DD→MN synaptic weight (µS)
# Optimization settings
N_TRIALS = 25 # Number of optimization trials (increase for production use)
TIMEOUT_SECONDS = 36000 # Maximum optimization time (10 hours)
# Gfluctdv (fluctuating conductance noise) parameters
ENABLE_GFLUCTDV = False # Enable membrane noise in motor neurons
GFLUCTDV_NOISE_MIN = 5e-6 # Minimum noise amplitude (S/cm²)
GFLUCTDV_NOISE_MAX = 3e-5 # Maximum noise amplitude (S/cm²)
# Results directory - use project root to share across examples
# Handle both manual execution and Sphinx Gallery (where __file__ is not defined)
try:
_script_dir = Path(__file__).parent
except NameError:
_script_dir = Path.cwd()
RESULTS_DIR = _script_dir.parent.parent / "results" / "dd_optimization"
RESULTS_DIR.mkdir(exist_ok=True, parents=True)
Load NEURON Mechanisms#
Compile and load the custom NMODL mechanisms needed for biophysical simulation.
load_nmodl_mechanisms()
h.secondorder = 2 # Use Crank-Nicolson integration (second-order accurate)
Define Objective Function#
The objective function evaluates network performance for a given parameter set. Optuna will minimize this function to find optimal parameters.
Optimization variables (suggested by Optuna each trial):
dd_neurons: Number of descending drive neurons [100-1000]conn_probability: DD→MN connection probability [0.1-1.0]dd_drive__Hz: Input drive frequency [5-1000 Hz]gfluctdv_noise_amplitude: Membrane noise level (if enabled)
Objective components:
mean_error: Normalized deviation from target mean firing ratestd_error: Normalized deviation from target firing rate stdplausibility_penalty: Soft constraints on biologically unrealistic parameters
def objective(trial):
"""
Optuna single-objective optimization function.
Optimizes network parameters to match target firing rate statistics
while preferring biologically plausible configurations.
Parameters
----------
trial : optuna.Trial
Current optimization trial with parameter suggestions
Returns
-------
float
Objective value (lower is better)
"""
try:
# Suggest parameters for this trial
dd_neurons = trial.suggest_int("dd_neurons", N_DD_NEURONS_MIN, N_DD_NEURONS_MAX)
conn_probability = trial.suggest_float("conn_prob", 0.1, 1.0)
dd_drive__Hz = trial.suggest_float("dd_drive", 5.0, 1000.0)
# Use fixed gamma shape (not optimized)
gamma_shape = get_gamma_shape_for_mvc(100, mvc_shape_value=GAMMA_SHAPE)
# Optional: Gfluctdv noise amplitude
if ENABLE_GFLUCTDV:
gfluctdv_noise_amplitude = trial.suggest_float(
"gfluctdv_noise_amplitude", GFLUCTDV_NOISE_MIN, GFLUCTDV_NOISE_MAX
)
else:
gfluctdv_noise_amplitude = 0.0
# Create recruitment thresholds for motor unit pool
recruitment_thresholds, _ = RecruitmentThresholds(
N=N_MOTOR_UNITS,
recruitment_range__ratio=100,
deluca__slope=5,
konstantin__max_threshold__ratio=1.0,
mode="combined",
)
# Create motor neuron pool
motor_neuron_pool = AlphaMN__Pool(
recruitment_thresholds__array=recruitment_thresholds,
config_file="alpha_mn_default.yaml",
)
# Apply Gfluctdv noise if enabled
if ENABLE_GFLUCTDV:
for cell in motor_neuron_pool:
cell.insert_Gfluctdv()
for d in cell.dend:
d.std_e_Gfluctdv = gfluctdv_noise_amplitude
d.std_i_Gfluctdv = gfluctdv_noise_amplitude
# Create descending drive pool
descending_drive_pool = DescendingDrive__Pool(
n=dd_neurons,
timestep__ms=TIMESTEP_MS * pq.ms,
process_type="gamma",
shape=gamma_shape, # type: ignore
)
# Build network with synaptic connections
network = Network({"DD": descending_drive_pool, "aMN": motor_neuron_pool})
network.connect(
source="DD",
target="aMN",
probability=conn_probability,
weight__uS=SYNAPTIC_WEIGHT * pq.uS,
)
network.connect_from_external(source="cortical_input", target="DD", weight__uS=1.0 * pq.uS)
dd_netcons = network.get_netcons("cortical_input", "DD")
# Setup spike recording
mn_spike_recorders = []
for cell in motor_neuron_pool:
spike_recorder = h.Vector()
nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
nc.threshold = 50
nc.record(spike_recorder)
mn_spike_recorders.append(spike_recorder)
# Generate constant drive signal with small noise
time_points = int(SIMULATION_TIME_MS / TIMESTEP_MS)
drive_signal = np.ones(time_points) * dd_drive__Hz + np.clip(
RANDOM_GENERATOR.normal(0, 1.0, size=time_points), 0, None
)
# Run NEURON simulation
h.load_file("stdrun.hoc")
h.dt = TIMESTEP_MS
h.tstop = SIMULATION_TIME_MS
# Initialize voltages
for section, voltage in zip(*motor_neuron_pool.get_initialization_data()):
section.v = voltage
for section, voltage in zip(*descending_drive_pool.get_initialization_data()):
section.v = voltage
h.finitialize()
# Simulation loop with DD drive injection
step_counter = 0
while h.t < h.tstop:
current_drive = drive_signal[min(step_counter, len(drive_signal) - 1)]
for dd_cell in descending_drive_pool:
if dd_cell.integrate(current_drive):
if h.t < h.tstop:
dd_netcons[dd_cell.pool__ID].event(h.t + 1)
h.fadvance()
step_counter += 1
# Convert spike data to Neo format
mn_segment = Segment(name="Motor Neurons")
dt_s = h.dt / 1000.0
mn_segment.spiketrains = [
SpikeTrain(
recorder.as_numpy() / 1000 * pq.s,
t_stop=SIMULATION_TIME_MS / 1000 * pq.s,
sampling_rate=(1 / dt_s * (pq.Hz)),
sampling_period=dt_s * pq.s,
name=f"MN_{i}",
)
for i, recorder in enumerate(mn_spike_recorders)
]
# Calculate firing rate statistics
stats = calculate_firing_rate_statistics(mn_segment.spiketrains)
n_active: int = int(stats["n_active"])
# Penalty for insufficient recruitment
if n_active < 10:
return 1000.0
fr_mean: float = float(stats["FR_mean"])
fr_std: float = float(stats["FR_std"])
# Normalized errors
mean_error = abs(fr_mean - TARGET_FR_MEAN__HZ) / TARGET_FR_MEAN__HZ
std_error = abs(fr_std - TARGET_FR_STD__HZ) / TARGET_FR_STD__HZ
# Biological plausibility penalties
plausibility_penalty = 0.0
if conn_probability > 0.7: # Very high connectivity is less realistic
plausibility_penalty += 0.1 * (conn_probability - 0.7)
if dd_neurons > 800: # Excessive neurons are less common
plausibility_penalty += 0.1 * ((dd_neurons - 800) / 200)
# Combined objective (minimize)
objective_value = mean_error + std_error + plausibility_penalty
# Store trial metadata
trial.set_user_attr("n_active", n_active)
trial.set_user_attr("FR_mean", fr_mean)
trial.set_user_attr("FR_std", fr_std)
trial.set_user_attr("mean_error", float(mean_error))
trial.set_user_attr("std_error", float(std_error))
trial.set_user_attr("plausibility_penalty", float(plausibility_penalty))
trial.set_user_attr("objective_value", float(objective_value))
trial.set_user_attr("dd_neurons", dd_neurons)
trial.set_user_attr("conn_probability", float(conn_probability))
trial.set_user_attr("synaptic_weight", float(SYNAPTIC_WEIGHT))
trial.set_user_attr("dd_drive__Hz", float(dd_drive__Hz))
trial.set_user_attr("gamma_shape", float(gamma_shape))
trial.set_user_attr("gfluctdv_enabled", ENABLE_GFLUCTDV)
trial.set_user_attr("gfluctdv_noise_amplitude", float(gfluctdv_noise_amplitude))
# Progress reporting
if trial.number % 1 == 0:
print(
f"Trial {trial.number}: FR={fr_mean:.1f}±{fr_std:.1f}Hz, "
f"obj={objective_value:.3f} (mean_err={mean_error:.3f}, std_err={std_error:.3f})"
)
return objective_value
except Exception as e:
print(f"Trial {trial.number} failed with error: {e}")
return 1000.0
Run Optimization Study#
Create an Optuna study and optimize parameters using Tree-structured Parzen Estimator (TPE). TPE is a Bayesian optimization algorithm that builds probabilistic models of the objective function.
print(
f"\nOptimizing network parameters | "
f"Target: {TARGET_FR_MEAN__HZ:.1f}±{TARGET_FR_STD__HZ:.1f}Hz | "
f"Trials: {N_TRIALS}\n"
)
# Create or load existing study (allows resuming interrupted optimization)
storage_name = f"sqlite:///{RESULTS_DIR}/optuna_dd_optimization.db"
study = optuna.create_study(
direction="minimize", # Minimize objective value
sampler=optuna.samplers.TPESampler(seed=42, multivariate=True),
study_name="dd_optimization",
storage=storage_name,
load_if_exists=True,
)
# Run optimization
study.optimize(objective, n_trials=N_TRIALS, timeout=TIMEOUT_SECONDS, show_progress_bar=True)
Optimizing network parameters | Target: 16.8±2.5Hz | Trials: 25
[I 2025-12-29 17:00:52,169] Using an existing study with name 'dd_optimization' instead of creating a new one.
0%| | 0/25 [00:00<?, ?it/s]Trial 25: FR=165.0±145.6Hz, obj=66.063 (mean_err=8.809, std_err=57.226)
[I 2025-12-29 17:01:01,360] Trial 25 finished with value: 66.06339822077369 and parameters: {'dd_neurons': 368, 'conn_prob': 0.9905094773673985, 'dd_drive': 387.8474799592962}. Best is trial 25 with value: 66.06339822077369.
0%| | 0/25 [00:09<?, ?it/s]
Best trial: 25. Best value: 66.0634: 0%| | 0/25 [00:09<?, ?it/s]
Best trial: 25. Best value: 66.0634: 4%|▍ | 1/25 [00:09<03:40, 9.19s/it]
Best trial: 25. Best value: 66.0634: 4%|▍ | 1/25 [00:09<03:40, 9.19s/it, 9.18/36000 seconds]Trial 26: FR=71.3±6.2Hz, obj=4.758 (mean_err=3.242, std_err=1.490)
[I 2025-12-29 17:01:06,529] Trial 26 finished with value: 4.758177967882577 and parameters: {'dd_neurons': 235, 'conn_prob': 0.9671907398425668, 'dd_drive': 259.81204655788366}. Best is trial 26 with value: 4.758177967882577.
Best trial: 25. Best value: 66.0634: 4%|▍ | 1/25 [00:14<03:40, 9.19s/it, 9.18/36000 seconds]
Best trial: 26. Best value: 4.75818: 4%|▍ | 1/25 [00:14<03:40, 9.19s/it, 9.18/36000 seconds]
Best trial: 26. Best value: 4.75818: 8%|▊ | 2/25 [00:14<02:36, 6.82s/it, 9.18/36000 seconds]
Best trial: 26. Best value: 4.75818: 8%|▊ | 2/25 [00:14<02:36, 6.82s/it, 14.35/36000 seconds]Trial 27: FR=19.4±3.3Hz, obj=0.487 (mean_err=0.154, std_err=0.310)
[I 2025-12-29 17:01:09,255] Trial 27 finished with value: 0.48655363546143326 and parameters: {'dd_neurons': 140, 'conn_prob': 0.9278382543576947, 'dd_drive': 83.35992660903196}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 26. Best value: 4.75818: 8%|▊ | 2/25 [00:17<02:36, 6.82s/it, 14.35/36000 seconds]
Best trial: 27. Best value: 0.486554: 8%|▊ | 2/25 [00:17<02:36, 6.82s/it, 14.35/36000 seconds]
Best trial: 27. Best value: 0.486554: 12%|█▏ | 3/25 [00:17<01:48, 4.95s/it, 14.35/36000 seconds]
Best trial: 27. Best value: 0.486554: 12%|█▏ | 3/25 [00:17<01:48, 4.95s/it, 17.08/36000 seconds]Trial 28: FR=81.7±6.8Hz, obj=5.591 (mean_err=3.855, std_err=1.708)
[I 2025-12-29 17:01:15,231] Trial 28 finished with value: 5.5914837384935465 and parameters: {'dd_neurons': 285, 'conn_prob': 0.9881725341563736, 'dd_drive': 255.47301752463798}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 12%|█▏ | 3/25 [00:23<01:48, 4.95s/it, 17.08/36000 seconds]
Best trial: 27. Best value: 0.486554: 12%|█▏ | 3/25 [00:23<01:48, 4.95s/it, 17.08/36000 seconds]
Best trial: 27. Best value: 0.486554: 16%|█▌ | 4/25 [00:23<01:52, 5.36s/it, 17.08/36000 seconds]
Best trial: 27. Best value: 0.486554: 16%|█▌ | 4/25 [00:23<01:52, 5.36s/it, 23.06/36000 seconds]Trial 29: FR=27.5±4.9Hz, obj=1.619 (mean_err=0.637, std_err=0.956)
[I 2025-12-29 17:01:18,659] Trial 29 finished with value: 1.6188587796479375 and parameters: {'dd_neurons': 211, 'conn_prob': 0.9522656887511345, 'dd_drive': 82.30407849542871}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 16%|█▌ | 4/25 [00:26<01:52, 5.36s/it, 23.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 16%|█▌ | 4/25 [00:26<01:52, 5.36s/it, 23.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 20%|██ | 5/25 [00:26<01:33, 4.66s/it, 23.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 20%|██ | 5/25 [00:26<01:33, 4.66s/it, 26.48/36000 seconds]Trial 30: FR=11.2±2.9Hz, obj=0.509 (mean_err=0.333, std_err=0.154)
[I 2025-12-29 17:01:22,237] Trial 30 finished with value: 0.5092425626551388 and parameters: {'dd_neurons': 291, 'conn_prob': 0.9131424928659169, 'dd_drive': 24.405849939824357}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 20%|██ | 5/25 [00:30<01:33, 4.66s/it, 26.48/36000 seconds]
Best trial: 27. Best value: 0.486554: 20%|██ | 5/25 [00:30<01:33, 4.66s/it, 26.48/36000 seconds]
Best trial: 27. Best value: 0.486554: 24%|██▍ | 6/25 [00:30<01:21, 4.29s/it, 26.48/36000 seconds]
Best trial: 27. Best value: 0.486554: 24%|██▍ | 6/25 [00:30<01:21, 4.29s/it, 30.06/36000 seconds]Trial 31: FR=6.4±1.9Hz, obj=0.876 (mean_err=0.620, std_err=0.251)
[I 2025-12-29 17:01:25,049] Trial 31 finished with value: 0.8758720442270053 and parameters: {'dd_neurons': 198, 'conn_prob': 0.7460283677095826, 'dd_drive': 32.25743637790005}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 24%|██▍ | 6/25 [00:32<01:21, 4.29s/it, 30.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 24%|██▍ | 6/25 [00:32<01:21, 4.29s/it, 30.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 28%|██▊ | 7/25 [00:32<01:08, 3.81s/it, 30.06/36000 seconds]
Best trial: 27. Best value: 0.486554: 28%|██▊ | 7/25 [00:32<01:08, 3.81s/it, 32.87/36000 seconds]Trial 32: FR=2.7±0.9Hz, obj=1.475 (mean_err=0.842, std_err=0.630)
[I 2025-12-29 17:01:28,305] Trial 32 finished with value: 1.4750543361207376 and parameters: {'dd_neurons': 271, 'conn_prob': 0.7299139787722873, 'dd_drive': 20.151877952780218}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 28%|██▊ | 7/25 [00:36<01:08, 3.81s/it, 32.87/36000 seconds]
Best trial: 27. Best value: 0.486554: 28%|██▊ | 7/25 [00:36<01:08, 3.81s/it, 32.87/36000 seconds]
Best trial: 27. Best value: 0.486554: 32%|███▏ | 8/25 [00:36<01:01, 3.63s/it, 32.87/36000 seconds]
Best trial: 27. Best value: 0.486554: 32%|███▏ | 8/25 [00:36<01:01, 3.63s/it, 36.13/36000 seconds]
[I 2025-12-29 17:01:30,537] Trial 33 finished with value: 1000.0 and parameters: {'dd_neurons': 141, 'conn_prob': 0.5911688331432965, 'dd_drive': 21.965126400246163}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 32%|███▏ | 8/25 [00:38<01:01, 3.63s/it, 36.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 32%|███▏ | 8/25 [00:38<01:01, 3.63s/it, 36.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 36%|███▌ | 9/25 [00:38<00:51, 3.19s/it, 36.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 36%|███▌ | 9/25 [00:38<00:51, 3.19s/it, 38.36/36000 seconds]Trial 34: FR=26.5±3.8Hz, obj=1.089 (mean_err=0.573, std_err=0.510)
[I 2025-12-29 17:01:34,135] Trial 34 finished with value: 1.0888004031118066 and parameters: {'dd_neurons': 245, 'conn_prob': 0.7547319579349063, 'dd_drive': 83.95874986642848}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 36%|███▌ | 9/25 [00:41<00:51, 3.19s/it, 38.36/36000 seconds]
Best trial: 27. Best value: 0.486554: 36%|███▌ | 9/25 [00:41<00:51, 3.19s/it, 38.36/36000 seconds]
Best trial: 27. Best value: 0.486554: 40%|████ | 10/25 [00:41<00:49, 3.32s/it, 38.36/36000 seconds]
Best trial: 27. Best value: 0.486554: 40%|████ | 10/25 [00:41<00:49, 3.32s/it, 41.96/36000 seconds]Trial 35: FR=8.2±2.2Hz, obj=0.661 (mean_err=0.514, std_err=0.122)
[I 2025-12-29 17:01:38,503] Trial 35 finished with value: 0.6612383732255603 and parameters: {'dd_neurons': 394, 'conn_prob': 0.9467755570026145, 'dd_drive': 13.995251228238617}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 40%|████ | 10/25 [00:46<00:49, 3.32s/it, 41.96/36000 seconds]
Best trial: 27. Best value: 0.486554: 40%|████ | 10/25 [00:46<00:49, 3.32s/it, 41.96/36000 seconds]
Best trial: 27. Best value: 0.486554: 44%|████▍ | 11/25 [00:46<00:50, 3.64s/it, 41.96/36000 seconds]
Best trial: 27. Best value: 0.486554: 44%|████▍ | 11/25 [00:46<00:50, 3.64s/it, 46.33/36000 seconds]Trial 36: FR=18.2±3.7Hz, obj=0.562 (mean_err=0.082, std_err=0.461)
[I 2025-12-29 17:01:43,248] Trial 36 finished with value: 0.5617209057627964 and parameters: {'dd_neurons': 442, 'conn_prob': 0.889277094300743, 'dd_drive': 25.830292228812947}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 44%|████▍ | 11/25 [00:51<00:50, 3.64s/it, 46.33/36000 seconds]
Best trial: 27. Best value: 0.486554: 44%|████▍ | 11/25 [00:51<00:50, 3.64s/it, 46.33/36000 seconds]
Best trial: 27. Best value: 0.486554: 48%|████▊ | 12/25 [00:51<00:51, 3.98s/it, 46.33/36000 seconds]
Best trial: 27. Best value: 0.486554: 48%|████▊ | 12/25 [00:51<00:51, 3.98s/it, 51.07/36000 seconds]
[I 2025-12-29 17:01:47,439] Trial 37 finished with value: 1000.0 and parameters: {'dd_neurons': 406, 'conn_prob': 0.9218021603221055, 'dd_drive': 10.088857037199556}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 48%|████▊ | 12/25 [00:55<00:51, 3.98s/it, 51.07/36000 seconds]
Best trial: 27. Best value: 0.486554: 48%|████▊ | 12/25 [00:55<00:51, 3.98s/it, 51.07/36000 seconds]
Best trial: 27. Best value: 0.486554: 52%|█████▏ | 13/25 [00:55<00:48, 4.04s/it, 51.07/36000 seconds]
Best trial: 27. Best value: 0.486554: 52%|█████▏ | 13/25 [00:55<00:48, 4.04s/it, 55.26/36000 seconds]Trial 38: FR=23.8±4.1Hz, obj=1.082 (mean_err=0.415, std_err=0.644)
[I 2025-12-29 17:01:53,493] Trial 38 finished with value: 1.0820509493050878 and parameters: {'dd_neurons': 611, 'conn_prob': 0.9336370371554951, 'dd_drive': 23.71658339107741}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 52%|█████▏ | 13/25 [01:01<00:48, 4.04s/it, 55.26/36000 seconds]
Best trial: 27. Best value: 0.486554: 52%|█████▏ | 13/25 [01:01<00:48, 4.04s/it, 55.26/36000 seconds]
Best trial: 27. Best value: 0.486554: 56%|█████▌ | 14/25 [01:01<00:51, 4.65s/it, 55.26/36000 seconds]
Best trial: 27. Best value: 0.486554: 56%|█████▌ | 14/25 [01:01<00:51, 4.65s/it, 61.32/36000 seconds]Trial 39: FR=37.7±5.7Hz, obj=2.549 (mean_err=1.242, std_err=1.278)
[I 2025-12-29 17:01:58,304] Trial 39 finished with value: 2.5491850364261155 and parameters: {'dd_neurons': 370, 'conn_prob': 0.9888315832434473, 'dd_drive': 67.14747905674382}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 56%|█████▌ | 14/25 [01:06<00:51, 4.65s/it, 61.32/36000 seconds]
Best trial: 27. Best value: 0.486554: 56%|█████▌ | 14/25 [01:06<00:51, 4.65s/it, 61.32/36000 seconds]
Best trial: 27. Best value: 0.486554: 60%|██████ | 15/25 [01:06<00:46, 4.70s/it, 61.32/36000 seconds]
Best trial: 27. Best value: 0.486554: 60%|██████ | 15/25 [01:06<00:46, 4.70s/it, 66.13/36000 seconds]
[I 2025-12-29 17:02:00,745] Trial 40 finished with value: 1000.0 and parameters: {'dd_neurons': 165, 'conn_prob': 0.9452930627694545, 'dd_drive': 13.035195361784247}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 60%|██████ | 15/25 [01:08<00:46, 4.70s/it, 66.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 60%|██████ | 15/25 [01:08<00:46, 4.70s/it, 66.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 64%|██████▍ | 16/25 [01:08<00:36, 4.02s/it, 66.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 64%|██████▍ | 16/25 [01:08<00:36, 4.02s/it, 68.57/36000 seconds]Trial 41: FR=21.4±4.1Hz, obj=0.900 (mean_err=0.272, std_err=0.622)
[I 2025-12-29 17:02:05,551] Trial 41 finished with value: 0.9001427102545785 and parameters: {'dd_neurons': 458, 'conn_prob': 0.7610005971749821, 'dd_drive': 34.93865690666536}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 64%|██████▍ | 16/25 [01:13<00:36, 4.02s/it, 68.57/36000 seconds]
Best trial: 27. Best value: 0.486554: 64%|██████▍ | 16/25 [01:13<00:36, 4.02s/it, 68.57/36000 seconds]
Best trial: 27. Best value: 0.486554: 68%|██████▊ | 17/25 [01:13<00:34, 4.26s/it, 68.57/36000 seconds]
Best trial: 27. Best value: 0.486554: 68%|██████▊ | 17/25 [01:13<00:34, 4.26s/it, 73.38/36000 seconds]Trial 42: FR=28.5±4.1Hz, obj=1.351 (mean_err=0.693, std_err=0.647)
[I 2025-12-29 17:02:08,573] Trial 42 finished with value: 1.3510958561388005 and parameters: {'dd_neurons': 161, 'conn_prob': 0.8035925873056585, 'dd_drive': 132.43934944638966}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 68%|██████▊ | 17/25 [01:16<00:34, 4.26s/it, 73.38/36000 seconds]
Best trial: 27. Best value: 0.486554: 68%|██████▊ | 17/25 [01:16<00:34, 4.26s/it, 73.38/36000 seconds]
Best trial: 27. Best value: 0.486554: 72%|███████▏ | 18/25 [01:16<00:27, 3.88s/it, 73.38/36000 seconds]
Best trial: 27. Best value: 0.486554: 72%|███████▏ | 18/25 [01:16<00:27, 3.88s/it, 76.40/36000 seconds]
[I 2025-12-29 17:02:12,612] Trial 43 finished with value: 1000.0 and parameters: {'dd_neurons': 409, 'conn_prob': 0.8073110075900845, 'dd_drive': 7.6513518696384715}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 72%|███████▏ | 18/25 [01:20<00:27, 3.88s/it, 76.40/36000 seconds]
Best trial: 27. Best value: 0.486554: 72%|███████▏ | 18/25 [01:20<00:27, 3.88s/it, 76.40/36000 seconds]
Best trial: 27. Best value: 0.486554: 76%|███████▌ | 19/25 [01:20<00:23, 3.93s/it, 76.40/36000 seconds]
Best trial: 27. Best value: 0.486554: 76%|███████▌ | 19/25 [01:20<00:23, 3.93s/it, 80.44/36000 seconds]Trial 44: FR=86.5±7.0Hz, obj=5.977 (mean_err=4.141, std_err=1.811)
[I 2025-12-29 17:02:20,310] Trial 44 finished with value: 5.976675501094426 and parameters: {'dd_neurons': 491, 'conn_prob': 0.9421420864108178, 'dd_drive': 169.5538515858688}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 76%|███████▌ | 19/25 [01:28<00:23, 3.93s/it, 80.44/36000 seconds]
Best trial: 27. Best value: 0.486554: 76%|███████▌ | 19/25 [01:28<00:23, 3.93s/it, 80.44/36000 seconds]
Best trial: 27. Best value: 0.486554: 80%|████████ | 20/25 [01:28<00:25, 5.06s/it, 80.44/36000 seconds]
Best trial: 27. Best value: 0.486554: 80%|████████ | 20/25 [01:28<00:25, 5.06s/it, 88.13/36000 seconds]
[I 2025-12-29 17:02:22,455] Trial 45 finished with value: 1000.0 and parameters: {'dd_neurons': 130, 'conn_prob': 0.8759934690623258, 'dd_drive': 13.001602761085465}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 80%|████████ | 20/25 [01:30<00:25, 5.06s/it, 88.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 80%|████████ | 20/25 [01:30<00:25, 5.06s/it, 88.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 84%|████████▍ | 21/25 [01:30<00:16, 4.19s/it, 88.13/36000 seconds]
Best trial: 27. Best value: 0.486554: 84%|████████▍ | 21/25 [01:30<00:16, 4.19s/it, 90.28/36000 seconds]Trial 46: FR=7.9±2.3Hz, obj=0.610 (mean_err=0.529, std_err=0.063)
[I 2025-12-29 17:02:26,337] Trial 46 finished with value: 0.6097764652957746 and parameters: {'dd_neurons': 356, 'conn_prob': 0.8760726858220766, 'dd_drive': 16.815866371005477}. Best is trial 27 with value: 0.48655363546143326.
Best trial: 27. Best value: 0.486554: 84%|████████▍ | 21/25 [01:34<00:16, 4.19s/it, 90.28/36000 seconds]
Best trial: 27. Best value: 0.486554: 84%|████████▍ | 21/25 [01:34<00:16, 4.19s/it, 90.28/36000 seconds]
Best trial: 27. Best value: 0.486554: 88%|████████▊ | 22/25 [01:34<00:12, 4.10s/it, 90.28/36000 seconds]
Best trial: 27. Best value: 0.486554: 88%|████████▊ | 22/25 [01:34<00:12, 4.10s/it, 94.16/36000 seconds]Trial 47: FR=10.9±2.4Hz, obj=0.414 (mean_err=0.355, std_err=0.045)
[I 2025-12-29 17:02:30,416] Trial 47 finished with value: 0.41350330058970947 and parameters: {'dd_neurons': 379, 'conn_prob': 0.841484095771859, 'dd_drive': 19.433516412986123}. Best is trial 47 with value: 0.41350330058970947.
Best trial: 27. Best value: 0.486554: 88%|████████▊ | 22/25 [01:38<00:12, 4.10s/it, 94.16/36000 seconds]
Best trial: 47. Best value: 0.413503: 88%|████████▊ | 22/25 [01:38<00:12, 4.10s/it, 94.16/36000 seconds]
Best trial: 47. Best value: 0.413503: 92%|█████████▏| 23/25 [01:38<00:08, 4.09s/it, 94.16/36000 seconds]
Best trial: 47. Best value: 0.413503: 92%|█████████▏| 23/25 [01:38<00:08, 4.09s/it, 98.24/36000 seconds]Trial 48: FR=51.3±6.3Hz, obj=3.592 (mean_err=2.048, std_err=1.529)
[I 2025-12-29 17:02:35,657] Trial 48 finished with value: 3.591537143051338 and parameters: {'dd_neurons': 368, 'conn_prob': 0.8474686543048925, 'dd_drive': 119.35463235652679}. Best is trial 47 with value: 0.41350330058970947.
Best trial: 47. Best value: 0.413503: 92%|█████████▏| 23/25 [01:43<00:08, 4.09s/it, 98.24/36000 seconds]
Best trial: 47. Best value: 0.413503: 92%|█████████▏| 23/25 [01:43<00:08, 4.09s/it, 98.24/36000 seconds]
Best trial: 47. Best value: 0.413503: 96%|█████████▌| 24/25 [01:43<00:04, 4.44s/it, 98.24/36000 seconds]
Best trial: 47. Best value: 0.413503: 96%|█████████▌| 24/25 [01:43<00:04, 4.44s/it, 103.48/36000 seconds]Trial 49: FR=10.3±3.7Hz, obj=0.859 (mean_err=0.387, std_err=0.473)
[I 2025-12-29 17:02:38,559] Trial 49 finished with value: 0.8594808661111732 and parameters: {'dd_neurons': 124, 'conn_prob': 0.10989027031939413, 'dd_drive': 435.5300626213635}. Best is trial 47 with value: 0.41350330058970947.
Best trial: 47. Best value: 0.413503: 96%|█████████▌| 24/25 [01:46<00:04, 4.44s/it, 103.48/36000 seconds]
Best trial: 47. Best value: 0.413503: 96%|█████████▌| 24/25 [01:46<00:04, 4.44s/it, 103.48/36000 seconds]
Best trial: 47. Best value: 0.413503: 100%|██████████| 25/25 [01:46<00:00, 3.98s/it, 103.48/36000 seconds]
Best trial: 47. Best value: 0.413503: 100%|██████████| 25/25 [01:46<00:00, 3.98s/it, 106.38/36000 seconds]
Best trial: 47. Best value: 0.413503: 100%|██████████| 25/25 [01:46<00:00, 4.26s/it, 106.38/36000 seconds]
Analyze Results#
Extract the best trial and display optimized parameters.
best_trial = study.best_trial
print(f"\nCompleted {len(study.trials)} trials")
print(
f"\nBest Trial {best_trial.number}: "
f"FR={best_trial.user_attrs.get('FR_mean'):.1f}±{best_trial.user_attrs.get('FR_std'):.1f}Hz | "
f"Objective={best_trial.value:.3f} | "
f"Drive={best_trial.user_attrs.get('dd_drive__Hz'):.1f}Hz | "
f"Neurons={best_trial.user_attrs.get('dd_neurons')} | "
f"Conn={best_trial.user_attrs.get('conn_probability'):.2f}"
)
Completed 50 trials
Best Trial 47: FR=10.9±2.4Hz | Objective=0.414 | Drive=19.4Hz | Neurons=379 | Conn=0.84
Save Optimized Parameters#
Store the best parameters in JSON format for use in production simulations.
def trial_to_dict(t):
"""Convert trial to dictionary of parameters."""
base_params = {
k: t.user_attrs.get(k)
for k in [
"FR_mean",
"FR_std",
"dd_neurons",
"conn_probability",
"dd_drive__Hz",
"gamma_shape",
"mean_error",
"std_error",
"plausibility_penalty",
"objective_value",
]
}
if ENABLE_GFLUCTDV:
base_params["gfluctdv_noise_amplitude"] = t.user_attrs.get("gfluctdv_noise_amplitude")
return base_params
results = {
"target": {
"FR_mean__Hz": TARGET_FR_MEAN__HZ,
"FR_std__Hz": TARGET_FR_STD__HZ,
},
"input_parameters": {
"gamma_shape": GAMMA_SHAPE,
"gfluctdv_enabled": ENABLE_GFLUCTDV,
},
"best_trial": trial_to_dict(best_trial),
}
json_path = RESULTS_DIR / "dd_optimized_params.json"
with open(json_path, "w") as f:
json.dump(results, f, indent=2)
joblib.dump(study, RESULTS_DIR / "study.pkl")
print(f"Saved optimized parameters: {json_path}\n")
Saved optimized parameters: /home/runner/work/MyoGen/MyoGen/results/dd_optimization/dd_optimized_params.json
Visualize Optimization History#
Plot the optimization trajectory showing how the objective value improved over trials.
# Get colors from the style
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
fig = plt.figure(figsize=(14, 10))
gs = fig.add_gridspec(2, 2, hspace=0.3, wspace=0.3)
# Left column: Optimization progress and FR convergence (shared x-axis)
ax_obj = fig.add_subplot(gs[0, 0])
ax_fr = fig.add_subplot(gs[1, 0], sharex=ax_obj)
# Right column: Parameter exploration and drive distribution
ax_param = fig.add_subplot(gs[0, 1])
ax_drive = fig.add_subplot(gs[1, 1])
# 1. Objective value over trials (log scale)
trial_numbers = [t.number for t in study.trials]
objective_values = [t.value for t in study.trials]
ax_obj.plot(trial_numbers, objective_values, "o-", alpha=0.6, markersize=4, color=colors[0])
ax_obj.axhline(best_trial.value, linestyle="--", color=colors[1], label=f"Best: {best_trial.value:.3f}")
ax_obj.set_ylabel("Objective Value (log scale)")
ax_obj.set_title("Optimization Progress")
ax_obj.set_yscale("log")
ax_obj.legend(framealpha=1.0, edgecolor="none")
ax_obj.tick_params(labelbottom=False)
# 2. Firing rate convergence
fr_means = [t.user_attrs.get("FR_mean", None) for t in study.trials]
fr_stds = [t.user_attrs.get("FR_std", None) for t in study.trials]
ax_fr.plot(trial_numbers, fr_means, "o-", alpha=0.6, markersize=4, color=colors[0], label="Mean FR")
ax_fr.axhline(TARGET_FR_MEAN__HZ, linestyle="--", color=colors[1], label=f"Target: {TARGET_FR_MEAN__HZ:.1f}Hz")
ax_fr.fill_between(
trial_numbers,
TARGET_FR_MEAN__HZ - TARGET_FR_STD__HZ,
TARGET_FR_MEAN__HZ + TARGET_FR_STD__HZ,
alpha=0.2,
color=colors[1],
label="Target Range",
)
ax_fr.set_xlabel("Trial")
ax_fr.set_ylabel("Firing Rate (Hz)")
ax_fr.set_title("Firing Rate Convergence")
ax_fr.legend(framealpha=1.0, edgecolor="none")
# 3. Parameter exploration: DD neurons vs Connection probability
dd_neurons_vals = [t.params.get("dd_neurons", None) for t in study.trials]
conn_prob_vals = [t.params.get("conn_prob", None) for t in study.trials]
scatter = ax_param.scatter(dd_neurons_vals, conn_prob_vals, c=objective_values, s=50, alpha=0.6)
ax_param.scatter(
best_trial.params["dd_neurons"],
best_trial.params["conn_prob"],
s=200,
marker="*",
color=colors[3],
label="Best Trial",
)
ax_param.set_xlabel("DD Neurons")
ax_param.set_ylabel("Connection Probability")
ax_param.set_title("Parameter Space Exploration")
ax_param.legend(framealpha=1.0, edgecolor="none")
plt.colorbar(scatter, ax=ax_param, label="Objective Value")
# 4. Drive frequency distribution
dd_drive_vals = [t.params.get("dd_drive", None) for t in study.trials]
ax_drive.hist(dd_drive_vals, bins=30, alpha=0.7, color=colors[2])
ax_drive.axvline(
best_trial.params["dd_drive"],
linestyle="--",
color=colors[1],
label=f'Best: {best_trial.params["dd_drive"]:.1f}Hz',
)
ax_drive.set_xlabel("DD Drive Frequency (Hz)")
ax_drive.set_ylabel("Number of Trials")
ax_drive.set_title("Drive Frequency Distribution")
ax_drive.legend(framealpha=1.0, edgecolor="none")
plt.tight_layout()
plt.savefig(RESULTS_DIR / "optimization_history.png", dpi=150, bbox_inches="tight")
plt.show()
Compare Target vs Achieved#
Validate that the optimization successfully matched the target firing rate statistics.
fig, ax = plt.subplots(1, 1, figsize=(10, 6))
# Bar plot comparing target vs achieved
categories = ["Mean FR (Hz)", "Std FR (Hz)"]
targets = [TARGET_FR_MEAN__HZ, TARGET_FR_STD__HZ]
achieved = [
best_trial.user_attrs.get("FR_mean"),
best_trial.user_attrs.get("FR_std"),
]
x = np.arange(len(categories))
width = 0.35
bars1 = ax.bar(x - width / 2, targets, width, label="Target")
bars2 = ax.bar(x + width / 2, achieved, width, label="Achieved")
# Add value labels on bars
for bars in [bars1, bars2]:
for bar in bars:
height = bar.get_height()
ax.text(
bar.get_x() + bar.get_width() / 2.0,
height,
f"{height:.2f}",
ha="center",
va="bottom",
)
# Calculate percent errors
mean_error = abs(achieved[0] - targets[0]) / targets[0] * 100
std_error = abs(achieved[1] - targets[1]) / targets[1] * 100
ax.set_ylabel("Value")
ax.set_title(
f"Firing Rate Optimization: Target vs Achieved\n"
f"(Mean error: {mean_error:.1f}%, Std error: {std_error:.1f}%)"
)
ax.set_xticks(x)
ax.set_xticklabels(categories)
ax.legend(framealpha=1.0, edgecolor="none")
plt.tight_layout()
plt.savefig(RESULTS_DIR / "target_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
Using Optimized Parameters#
Important
The optimized parameters can be loaded and used in production simulations:
import json
from pathlib import Path
# Load optimized parameters
with open("results/dd_optimization/dd_optimized_params.json") as f:
params = json.load(f)
# Extract best values
dd_neurons = int(params["best_trial"]["dd_neurons"])
conn_prob = params["best_trial"]["conn_probability"]
dd_drive_Hz = params["best_trial"]["dd_drive__Hz"]
# Use in your simulation
descending_drive_pool = DescendingDrive__Pool(n=dd_neurons, ...)
network.connect("DD", "aMN", probability=conn_prob, ...)
print("\n[DONE] Optimization complete!")
print(f"Best parameters saved to: {json_path}")
print(f"Full study saved to: {RESULTS_DIR / 'study.pkl'}")
[DONE] Optimization complete!
Best parameters saved to: /home/runner/work/MyoGen/MyoGen/results/dd_optimization/dd_optimized_params.json
Full study saved to: /home/runner/work/MyoGen/MyoGen/results/dd_optimization/study.pkl
Total running time of the script: (1 minutes 48.094 seconds)