Note

Go to the end to download the full example code.

Compute Muscle Force from Spinal Network Spike Trains#

This example demonstrates muscle force generation from the Watanabe spinal network simulation. It loads pre-computed spike trains and converts them to realistic muscle force output using the Fuglevand force model with batched processing for memory efficiency.

Note

Force computation pipeline:

  1. Load spike train data from previous analysis (NEO Block format)

  2. Extract alpha motor neuron spike trains

  3. Setup Fuglevand force model with physiological parameters

  4. Process motor units in batches (150 units at a time) for memory efficiency

  5. Accumulate force contributions from all motor units

  6. Save force output as NEO Block with comprehensive metadata

Important

Prerequisites: This example requires outputs from:

  • 03_10pct_mvc_simulation.py: Simulation parameters

  • 04_load_and_analyze_results.py: NEO-formatted spike train data (watanabe_results_neo.pkl)

Key Features:

  • Batched processing: Handles 800 motor units efficiently by processing 150 at a time

  • Fuglevand force model: Physiologically realistic twitch dynamics with 100:1 force range

  • High sampling rate: 2048 Hz force output for detailed temporal analysis

  • Motor unit diversity: 2x contraction time range with 90 ms slowest twitch rise time

  • NEO format output: Complete metadata preservation for reproducible analysis

  • Memory efficient: Batch processing with explicit garbage collection

MyoGen Components Used:

  • ForceModel: The Fuglevand force model converts spike trains to muscle force. Key parameters:

    • recruitment_thresholds: Determines each motor unit’s twitch amplitude (larger threshold → larger force, following the “size principle”)

    • recording_frequency__Hz: Output sampling rate (2048 Hz here)

    • longest_duration_rise_time__ms: Slowest motor unit’s twitch rise time

    • contraction_time_range_factor: Ratio of slowest/fastest contraction times

    The model generates a normalized twitch waveform for each motor unit, then convolves it with the spike train to produce force contributions.

  • RecruitmentThresholds: Same thresholds used during simulation - ensures force model matches network.

How Force Generation Works:

  1. Each motor unit has a characteristic twitch: F(t) = P * (t/T) * exp(1 - t/T) where P = peak force (from recruitment threshold) and T = contraction time.

  2. Smaller motor units (low threshold): small P, short T → small, fast twitches Larger motor units (high threshold): large P, long T → large, slow twitches

  3. Total muscle force = sum of all motor unit forces (unfused tetanus at ~10-30 Hz)

Use Case: Generate realistic muscle force output from motor neuron spike trains for biomechanical analysis, force-EMG relationship studies, and validation against experimental force recordings from Watanabe and Kohn (2015).

Workflow Position: Step 5 of 6

Next Step: Run 06_visualize.py to generate publication-quality figures.

Import Libraries#

from pathlib import Path

import joblib
import neo
import numpy as np
import quantities as pq

from myogen import simulator
from myogen.simulator.core.force.force_model import ForceModel

Setup Paths and Load Parameters#

try:
    _script_dir = Path(__file__).parent
except NameError:
    _script_dir = Path.cwd()

save_path = _script_dir / "results"
save_path.mkdir(exist_ok=True)

spinal_results_path = save_path / "watanabe_results_neo.pkl"

# Load simulation parameters
params_file = save_path / "watanabe_simulation_params.pkl"

if not params_file.exists():
    raise FileNotFoundError(
        f"Simulation parameters file not found: {params_file}\n"
        "Please run 03_10pct_mvc_simulation.py first to generate the simulation data."
    )

sim_params = joblib.load(params_file)
dt = sim_params["dt__ms"]  # ms - Integration timestep
tstop = sim_params["tstop__ms"]  # ms - Total simulation duration
naMN = sim_params["naMN"]  # Number of alpha motor neurons
segment_duration__s = sim_params["segment_duration__s"]

print("=" * 70)
print("SIMULATION PARAMETERS")
print("=" * 70)
print(f"Segment duration: {segment_duration__s} s")
print(f"Total duration: {tstop / 1000} s ({tstop} ms)")
print(f"Timestep: {dt} ms")
print(f"Motor neurons: {naMN}")
print()

# Memory optimization: process motor units in batches
motor_unit_batch_size = 150

======================================================================
SIMULATION PARAMETERS
======================================================================
Segment duration: 5 s
Total duration: 15.0 s (15000.0 ms)
Timestep: 0.025 ms
Motor neurons: 800

Load Spike Train Results#

print("=" * 70)
print("LOADING SPIKE TRAIN DATA")
print("=" * 70)

with open(spinal_results_path, "rb") as f:
    results: neo.Block = joblib.load(f)

# Extract alpha motor neuron spike trains
aMN_results: neo.Segment = results.filter(name="aMN", container=True)[0]
aMN_spikes = aMN_results.spiketrains

# Get actual number of motor neurons with spike trains
n_actual_motor_neurons = len(aMN_spikes)
print(f"Loaded {n_actual_motor_neurons} motor neuron spike trains (out of {naMN} expected)")

======================================================================
LOADING SPIKE TRAIN DATA
======================================================================
Loaded 671 motor neuron spike trains (out of 800 expected)

Setup Force Model Parameters#

The Fuglevand force model generates realistic motor unit twitch dynamics

print("\n" + "=" * 70)
print("CONFIGURING FORCE MODEL")
print("=" * 70)

# Generate recruitment thresholds for motor neuron pool
# These determine motor unit properties (twitch amplitude, contraction time)
recruitment_thresholds, _ = simulator.RecruitmentThresholds(
    N=naMN,
    recruitment_range__ratio=100,
    mode="combined",
    deluca__slope=10,
)

force_params = {
    "recording_frequency__Hz": 2048 * pq.Hz,  # Sampling frequency for force output
    "longest_duration_rise_time__ms": 90.0 * pq.ms,  # Slowest motor unit twitch rise time
    "contraction_time_range_factor": 2,  # Spread of contraction times across motor units
}

print(f"Recording frequency: {force_params['recording_frequency__Hz']}")
print(f"Slowest twitch rise time: {force_params['longest_duration_rise_time__ms']}")
print(f"Contraction time range factor: {force_params['contraction_time_range_factor']}")
print(f"Recruitment range ratio: 100:1 (smallest to largest motor unit)")

======================================================================
CONFIGURING FORCE MODEL
======================================================================
Recording frequency: 2048.0 Hz
Slowest twitch rise time: 90.0 ms
Contraction time range factor: 2
Recruitment range ratio: 100:1 (smallest to largest motor unit)

Generate Force Output in Batches#

Process motor units in batches for memory efficiency with large populations

print("\n" + "=" * 70)
print("GENERATING FORCE OUTPUT (BATCHED PROCESSING)")
print("=" * 70)
print(f"Processing {n_actual_motor_neurons} motor units in batches of {motor_unit_batch_size}")

# Calculate number of batches based on actual motor neurons with spikes
n_batches = int(np.ceil(n_actual_motor_neurons / motor_unit_batch_size))

# Initialize force accumulator
recording_freq_Hz = force_params["recording_frequency__Hz"].rescale("Hz").magnitude
force_duration_samples = int(tstop / 1000 * recording_freq_Hz)
force_total = np.zeros(force_duration_samples)

print(f"Total batches: {n_batches}")
print(f"Force signal samples: {force_duration_samples}")
print(f"Force signal duration: {tstop / 1000:.1f} s")

for batch_idx in range(n_batches):
    # Define motor unit batch indices (based on actual spike trains available)
    mu_start = batch_idx * motor_unit_batch_size
    mu_end = min((batch_idx + 1) * motor_unit_batch_size, n_actual_motor_neurons)

    print(f"\nProcessing batch {batch_idx + 1}/{n_batches}:")
    print(f"  Motor units: {mu_start + 1}-{mu_end} ({mu_end - mu_start} units)")

    # Create batch-specific Block with subset of motor units
    batch_block = neo.Block(name=f"MU_Batch_{batch_idx}")
    batch_segment = neo.Segment(name=f"aMN_batch_{batch_idx}")

    # Add only spike trains for this batch of motor units
    for mu_idx in range(mu_start, mu_end):
        batch_segment.spiketrains.append(aMN_spikes[mu_idx])

    batch_block.segments.append(batch_segment)

    # Create force model for this batch
    batch_recruitment_thresholds = recruitment_thresholds[mu_start:mu_end]
    batch_force_model = ForceModel(
        recruitment_thresholds=batch_recruitment_thresholds, **force_params
    )

    # Generate force for this batch
    force_batch = batch_force_model.generate_force(spike_train__Block=batch_block)
    force_batch_array = force_batch.magnitude.squeeze()

    print(f"  Batch force shape: {force_batch_array.shape}")
    print(f"  Batch force range: [{force_batch_array.min():.3f}, {force_batch_array.max():.3f}]")

    # Accumulate force (sum contributions from all motor units)
    force_total += force_batch_array

    # Clean up to free memory
    del batch_block, batch_segment, batch_force_model, force_batch, force_batch_array

# Create final AnalogSignal
force_output = neo.AnalogSignal(
    force_total * pq.dimensionless,
    t_start=0 * pq.ms,
    sampling_rate=force_params["recording_frequency__Hz"],
)

print("\n" + "=" * 70)
print("FORCE GENERATION COMPLETE")
print("=" * 70)
print(f"Final force signal shape: {force_output.shape}")
print(f"Final force range: [{force_total.min():.3f}, {force_total.max():.3f}]")
print(f"Sampling rate: {force_output.sampling_rate}")
print(f"Duration: {float(force_output.duration.rescale('s')):.2f} s")

======================================================================
GENERATING FORCE OUTPUT (BATCHED PROCESSING)
======================================================================
Processing 671 motor units in batches of 150
Total batches: 5
Force signal samples: 30720
Force signal duration: 15.0 s

Processing batch 1/5:
  Motor units: 1-150 (150 units)

Pool 1 Twitch trains (Numba):   0%|          | 0/150 [00:00<?, ?MU/s]
Pool 1 Twitch trains (Numba):   1%|          | 1/150 [00:00<00:17,  8.31MU/s]
Pool 1 Twitch trains (Numba):  48%|████▊     | 72/150 [00:00<00:00, 387.41MU/s]
Pool 1 Twitch trains (Numba): 100%|██████████| 150/150 [00:00<00:00, 480.87MU/s]
  Batch force shape: (30720,)
  Batch force range: [0.000, 28.265]

Processing batch 2/5:
  Motor units: 151-300 (150 units)

Pool 1 Twitch trains (Numba):   0%|          | 0/150 [00:00<?, ?MU/s]
Pool 1 Twitch trains (Numba):   1%|          | 1/150 [00:00<00:17,  8.47MU/s]
Pool 1 Twitch trains (Numba):  47%|████▋     | 70/150 [00:00<00:00, 381.27MU/s]
Pool 1 Twitch trains (Numba): 100%|██████████| 150/150 [00:00<00:00, 482.06MU/s]
  Batch force shape: (30720,)
  Batch force range: [0.000, 21.931]

Processing batch 3/5:
  Motor units: 301-450 (150 units)

Pool 1 Twitch trains (Numba):   0%|          | 0/150 [00:00<?, ?MU/s]
Pool 1 Twitch trains (Numba):   1%|          | 1/150 [00:00<00:17,  8.51MU/s]
Pool 1 Twitch trains (Numba):  49%|████▊     | 73/150 [00:00<00:00, 396.43MU/s]
Pool 1 Twitch trains (Numba): 100%|██████████| 150/150 [00:00<00:00, 490.95MU/s]
  Batch force shape: (30720,)
  Batch force range: [0.000, 20.042]

Processing batch 4/5:
  Motor units: 451-600 (150 units)

Pool 1 Twitch trains (Numba):   0%|          | 0/150 [00:00<?, ?MU/s]
Pool 1 Twitch trains (Numba):   1%|          | 1/150 [00:00<00:17,  8.54MU/s]
Pool 1 Twitch trains (Numba):  54%|█████▍    | 81/150 [00:00<00:00, 441.80MU/s]
Pool 1 Twitch trains (Numba): 100%|██████████| 150/150 [00:00<00:00, 543.31MU/s]
  Batch force shape: (30720,)
  Batch force range: [0.000, 17.431]

Processing batch 5/5:
  Motor units: 601-671 (71 units)

Pool 1 Twitch trains (Numba):   0%|          | 0/71 [00:00<?, ?MU/s]
Pool 1 Twitch trains (Numba):   1%|▏         | 1/71 [00:00<00:08,  8.58MU/s]
Pool 1 Twitch trains (Numba): 100%|██████████| 71/71 [00:00<00:00, 454.21MU/s]
  Batch force shape: (30720,)
  Batch force range: [0.000, 2.262]

======================================================================
FORCE GENERATION COMPLETE
======================================================================
Final force signal shape: (30720, 1)
Final force range: [0.000, 89.421]
Sampling rate: 2048.0 Hz
Duration: 15.00 s

Save Force Results as NEO Block#

Store force output with comprehensive metadata for reproducible analysis

print("\n" + "=" * 70)
print("SAVING FORCE RESULTS")
print("=" * 70)

# Create Neo Block to store force results with metadata
force_block = neo.Block(name="Watanabe Force Results")
force_segment = neo.Segment(name="Force")

# Add force output (already a Neo AnalogSignal)
force_segment.analogsignals.append(force_output)

# Store metadata as annotations
force_block.annotations["simulation_params"] = {
    "dt_ms": dt,
    "tstop_ms": tstop,
    "n_motor_neurons": naMN,
    "segment_duration_s": segment_duration__s,
}

force_block.annotations["recruitment_thresholds"] = recruitment_thresholds
force_block.annotations["force_model_params"] = force_params

# Create a reference force model for stats (using only motor units with spikes)
reference_force_model = ForceModel(
    recruitment_thresholds=recruitment_thresholds[:n_actual_motor_neurons],
    **force_params,
)

force_block.annotations["force_model_stats"] = {
    "n_motor_units": reference_force_model._number_of_neurons,
    "n_motor_units_with_spikes": n_actual_motor_neurons,
    "recruitment_ratio": float(reference_force_model._recruitment_ratio),
    "peak_force_range": [
        float(reference_force_model.peak_twitch_forces__unitless[0]),
        float(reference_force_model.peak_twitch_forces__unitless[-1]),
    ],
    "contraction_time_range_samples": [
        float(reference_force_model.contraction_times__samples[0]),
        float(reference_force_model.contraction_times__samples[-1]),
    ],
}

# Add segment to block
force_block.segments.append(force_segment)

# Save using joblib (compatible with Neo)
output_file = save_path / "watanabe__force_results.pkl"

with open(output_file, "wb") as f:
    joblib.dump(force_block, f)

print(f"Force results saved to: {output_file}")
print(f"File size: {output_file.stat().st_size / 1e6:.2f} MB")

======================================================================
SAVING FORCE RESULTS
======================================================================
Force results saved to: /home/runner/work/MyoGen/MyoGen/examples/03_papers/watanabe/results/watanabe__force_results.pkl
File size: 0.25 MB

Total running time of the script: (0 minutes 19.967 seconds)

Gallery generated by Sphinx-Gallery