Skip to content

Neuron injection & I/O

Current injection

inject_currents_into_populations 

inject_currents_into_populations(populations: Sequence[_Pool], input_current__AnalogSignal: CURRENT__AnalogSignal) -> None

Injects input currents into the specified populations.

Sets up time-varying current injection using NEURON's IClamp and Vector.play() mechanisms. This function only sets up the current injection - spike recording and simulation execution must be handled separately by the user.

Parameters:

Name Type Description Default
populations Sequence[_Pool]

The populations of neurons to inject current into.

 required
input_current__AnalogSignal CURRENT__AnalogSignal

The analog signal of input currents to inject into the population. Shape should be (time_points, n_pools) where n_pools matches len(populations).

 required

Returns:

Type Description
None

Current injection mechanisms are attached to the neurons as side effects.

Raises:

Type Description
ValueError

If the number of populations does not match the number of current channels.
Notes
  • Current injection vectors and IClamp objects are stored on each cell as cell._stim_vectors to prevent garbage collection
  • The user is responsible for setting up spike recording and calling h.run()
Source code in myogen/utils/neuron/inject_currents_into_populations.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
@beartowertype
def inject_currents_into_populations(
    populations: Sequence[_Pool],
    input_current__AnalogSignal: CURRENT__AnalogSignal,
) -> None:
    """
    Injects input currents into the specified populations.

    Sets up time-varying current injection using NEURON's IClamp and Vector.play()
    mechanisms. This function only sets up the current injection - spike recording
    and simulation execution must be handled separately by the user.

    Parameters
    ----------
    populations : Sequence[_Pool]
        The populations of neurons to inject current into.
    input_current__AnalogSignal : CURRENT__AnalogSignal
        The analog signal of input currents to inject into the population.
        Shape should be (time_points, n_pools) where n_pools matches len(populations).

    Returns
    -------
    None
        Current injection mechanisms are attached to the neurons as side effects.

    Raises
    ------
    ValueError
        If the number of populations does not match the number of current channels.

    Notes
    -----
    - Current injection vectors and IClamp objects are stored on each cell as
      `cell._stim_vectors` to prevent garbage collection
    - The user is responsible for setting up spike recording and calling h.run()
    """
    # Validate input dimensions
    n_pools = len(populations)
    n_current_channels = input_current__AnalogSignal.shape[1]

    if n_pools != n_current_channels:
        raise ValueError(
            f"Number of populations ({n_pools}) must match number of current channels "
            f"({n_current_channels})"
        )

    simulation_time__ms = input_current__AnalogSignal.t_stop.rescale(pq.ms).magnitude

    for pool_idx, pool in enumerate(populations):
        for cell in pool:
            # Create IClamp with time-varying amplitude using Vector.play()
            stim = h.IClamp(cell.soma(0.5))
            stim.delay = 0
            stim.dur = simulation_time__ms  # Keep it "always on"
            stim.amp = 0  # Will be controlled by vector

            # Create time and current vectors
            vec_t = h.Vector(input_current__AnalogSignal.times.rescale(pq.ms).magnitude)
            vec_i = h.Vector(
                input_current__AnalogSignal[:, pool_idx].rescale(pq.nA).magnitude
            )

            # Play waveform into stim.amp
            vec_i.play(stim._ref_amp, vec_t, 1)  # 1 = continuous update

            # CRITICAL: Keep references to vectors to prevent garbage collection
            # Store them as attributes on the cell to keep them alive
            cell._stim_vectors = (stim, vec_t, vec_i)  # Prevent GC cleanup

inject_currents_and_simulate_spike_trains 

inject_currents_and_simulate_spike_trains(populations: Sequence[_Pool], input_current__AnalogSignal: CURRENT__AnalogSignal, spike_detection_thresholds__mV: Quantity__mV | Sequence[Quantity__mV] = -10.0 * mV) -> SPIKE_TRAIN__Block

Injects input currents into populations and returns recorded spike trains.

This is a complete pipeline function that sets up current injection, spike recording, runs the NEURON simulation, and returns the results as a properly formatted neo.Block.

Parameters:

Name Type Description Default
populations Sequence[_Pool]

The populations of neurons to inject current into.

 required
input_current__AnalogSignal CURRENT__AnalogSignal

The analog signal of input currents to inject into the population. Shape should be (time_points, n_pools) where n_pools matches len(populations).

 required
spike_detection_thresholds__mV float | Sequence[float]

Thresholds for spike detection in millivolts, by default -10.0. If a sequence is provided, it must match the number of populations.

 -10.0 * mV

Returns:

Type Description
SPIKE_TRAIN__Block

Neo Block containing spike trains organized as segments (pools) with spiketrains (neurons). Each segment represents a motor unit pool, each spiketrain represents a neuron.

Raises:

Type Description
ValueError

If the number of populations does not match the number of current channels.
Notes

This function performs the complete simulation pipeline: 1. Sets up current injection (same as inject_currents_into_populations) 2. Sets up spike recording for all neurons 3. Runs the NEURON simulation via h.run() 4. Converts recorded spikes to neo.Block format

Source code in myogen/utils/neuron/inject_currents_into_populations.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
@beartowertype
def inject_currents_and_simulate_spike_trains(
    populations: Sequence[_Pool],
    input_current__AnalogSignal: CURRENT__AnalogSignal,
    spike_detection_thresholds__mV: Quantity__mV | Sequence[Quantity__mV] = -10.0
    * pq.mV,
) -> SPIKE_TRAIN__Block:
    """
    Injects input currents into populations and returns recorded spike trains.

    This is a complete pipeline function that sets up current injection, spike recording,
    runs the NEURON simulation, and returns the results as a properly formatted neo.Block.

    Parameters
    ----------
    populations : Sequence[_Pool]
        The populations of neurons to inject current into.
    input_current__AnalogSignal : CURRENT__AnalogSignal
        The analog signal of input currents to inject into the population.
        Shape should be (time_points, n_pools) where n_pools matches len(populations).
    spike_detection_thresholds__mV : float | Sequence[float], optional
        Thresholds for spike detection in millivolts, by default -10.0. If a sequence is provided, it must match the number of populations.


    Returns
    -------
    SPIKE_TRAIN__Block
        Neo Block containing spike trains organized as segments (pools) with spiketrains (neurons).
        Each segment represents a motor unit pool, each spiketrain represents a neuron.

    Raises
    ------
    ValueError
        If the number of populations does not match the number of current channels.

    Notes
    -----
    This function performs the complete simulation pipeline:
    1. Sets up current injection (same as inject_currents_into_populations)
    2. Sets up spike recording for all neurons
    3. Runs the NEURON simulation via h.run()
    4. Converts recorded spikes to neo.Block format
    """
    from neo import Block, Segment, SpikeTrain

    # First, set up current injection using the existing function
    inject_currents_into_populations(populations, input_current__AnalogSignal)

    simulation_time__ms = input_current__AnalogSignal.t_stop.rescale(pq.ms)

    # Validate input dimensions
    n_pools = len(populations)
    n_current_channels = input_current__AnalogSignal.shape[1]

    if n_pools != n_current_channels:
        raise ValueError(
            f"Number of populations ({n_pools}) must match number of current channels "
            f"({n_current_channels})"
        )

    if not isinstance(spike_detection_thresholds__mV, Sequence):
        spike_detection_thresholds__mV = [spike_detection_thresholds__mV] * n_pools

    # Set up spike recording for all neurons
    spike_recorders = []

    for pool_idx, pool in enumerate(populations):
        pool_spike_recorders = []

        for cell in pool:
            # Setup spike recording
            spike_recorder = h.Vector()
            nc = h.NetCon(cell.soma(0.5)._ref_v, None, sec=cell.soma)
            nc.threshold = spike_detection_thresholds__mV[pool_idx]
            nc.record(spike_recorder)

            pool_spike_recorders.append(spike_recorder)

        spike_recorders.append(pool_spike_recorders)

    # Initialize sections with proper voltages for each population
    for pool in populations:
        for section, voltage in zip(*pool.get_initialization_data()):
            section.v = voltage

    # Initialize and run the NEURON simulation
    h.finitialize()  # Use default initialization, voltages already set above
    # Advance the solver to the end of the simulation (replaces the deprecated
    # neuron.run(); finitialize() above already set the initial state).
    h.tstop = float(simulation_time__ms)
    while h.t < h.tstop:
        h.fadvance()

    # Convert spike data to neo.Block format
    block = Block()

    for pool_idx, pool_spike_recorders in enumerate(spike_recorders):
        # Create segment for this motor unit pool
        segment = Segment(name=f"Pool {pool_idx}")

        segment.spiketrains = [
            SpikeTrain(
                (spike_recorder.as_numpy() * pq.ms).rescale(pq.s),
                t_stop=simulation_time__ms.rescale(pq.s),
                sampling_rate=(1 / (h.dt * pq.ms)).rescale(pq.Hz),
                sampling_period=(h.dt * pq.ms).rescale(pq.s),
                name=str(int(neuron_idx)),
                description=f"Pool {pool_idx}, Neuron {neuron_idx}",
            )
            for neuron_idx, spike_recorder in enumerate(pool_spike_recorders)
        ]

        # Only add segment if it has spiketrains (which it always will now)
        block.segments.append(segment)

    return block

Persistence

ContinuousSaver 

ContinuousSaver(save_path: Path, chunk_duration__ms: Quantity__ms = 10000.0 * ms, populations: Optional[dict] = None, recording_config: Optional[dict] = None, verbose: bool = True)

Manages continuous saving of simulation data in chunks to prevent memory overflow.

Instead of accumulating all data in RAM, this class periodically saves chunks to disk and clears memory. Data can be loaded and combined afterward.

Parameters:

Name Type Description Default
save_path Path

Directory where chunks will be saved

 required
chunk_duration__ms float

Duration of each chunk in milliseconds (default: 10000 ms = 10 seconds)

 10000.0 * ms
populations dict

Dictionary of populations to record from

 None
recording_config dict

Configuration like {"aMN": [0, 10, 20, ...]} for which cells to record

 None
Source code in myogen/utils/continuous_saver.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def __init__(
    self,
    save_path: Path,
    chunk_duration__ms: Quantity__ms = 10000.0 * pq.ms,
    populations: Optional[dict] = None,
    recording_config: Optional[dict] = None,
    verbose: bool = True,
):
    self.save_path = Path(save_path)
    self.save_path.mkdir(exist_ok=True, parents=True)
    self.chunk_duration__ms = chunk_duration__ms
    self.populations = populations or {}
    self.recording_config = recording_config or {}
    self.verbose = verbose

    # Tracking state
    self.chunk_id = 0
    self.last_save_time = 0.0
    self.current_chunk_data = defaultdict(lambda: defaultdict(list))
    self.current_chunk_times = []

    # Spike recording
    self.spike_data = defaultdict(lambda: {"times": [], "ids": []})

    if self.verbose:
        print("ContinuousSaver initialized:")
        print(f"\tSave path: {self.save_path}")
        print(f"\tChunk duration: {chunk_duration__ms} ms")
        print(f"\tRecording config: {recording_config}")

record_step 

record_step(timestep__ms: float) -> None

Record data for current simulation timestep.

Call this from your step callback at each timestep.

Parameters:

Name Type Description Default
timestep__ms float

Integration timestep in milliseconds

 required
Source code in myogen/utils/continuous_saver.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def record_step(self, timestep__ms: float) -> None:
    """
    Record data for current simulation timestep.

    Call this from your step callback at each timestep.

    Parameters
    ----------
    timestep__ms : float
        Integration timestep in milliseconds
    """
    current_time = h.t

    # Record current time
    self.current_chunk_times.append(current_time)

    # Record membrane potentials for specified cells
    for pop_name, cell_indices in self.recording_config.items():
        if pop_name not in self.populations:
            continue

        population = self.populations[pop_name]

        for cell_idx in cell_indices:
            if cell_idx < len(population):
                # Read voltage from soma
                voltage = population[cell_idx].soma(0.5).v
                self.current_chunk_data[pop_name][cell_idx].append(voltage)

    # Check if it's time to save this chunk
    if current_time - self.last_save_time >= self.chunk_duration__ms:
        self._save_current_chunk(timestep__ms)

record_spike 

record_spike(pop_name: str, cell_id: int, spike_time: float) -> None

Record a spike event.

Parameters:

Name Type Description Default
pop_name str

Population name

 required
cell_id int

Cell ID within population

 required
spike_time float

Time of spike in milliseconds

 required
Source code in myogen/utils/continuous_saver.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def record_spike(self, pop_name: str, cell_id: int, spike_time: float) -> None:
    """
    Record a spike event.

    Parameters
    ----------
    pop_name : str
        Population name
    cell_id : int
        Cell ID within population
    spike_time : float
        Time of spike in milliseconds
    """
    self.spike_data[pop_name]["times"].append(spike_time)
    self.spike_data[pop_name]["ids"].append(cell_id)

finalize 

finalize(timestep__ms: Quantity__ms, spike_results=None) -> None

Save final chunk and spike data.

Call this after simulation completes.

Parameters:

Name Type Description Default
timestep__ms Quantity__ms

Integration timestep in milliseconds

 required
spike_results NEO Block

NEO Block containing spike trains from SimulationRunner. If provided, spike data will be extracted and saved to chunks.

 None
Source code in myogen/utils/continuous_saver.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
def finalize(self, timestep__ms: Quantity__ms, spike_results=None) -> None:
    """
    Save final chunk and spike data.

    Call this after simulation completes.

    Parameters
    ----------
    timestep__ms : Quantity__ms
        Integration timestep in milliseconds
    spike_results : NEO Block, optional
        NEO Block containing spike trains from SimulationRunner.
        If provided, spike data will be extracted and saved to chunks.
    """
    # Save any remaining data
    if len(self.current_chunk_times) > 0:
        self._save_current_chunk(timestep__ms)

    # Save spike data
    spike_filename = self.save_path / "spikes.pkl"
    spike_data_arrays = {}

    if spike_results is not None:
        # Extract spike data from NEO Block (from SimulationRunner)
        if self.verbose:
            print("\nExtracting spike data from SimulationRunner results...")
        from neo import Block

        if isinstance(spike_results, Block):
            for seg in spike_results.segments:
                if len(seg.spiketrains) > 0:
                    pop_name = seg.name
                    times_list = []
                    ids_list = []

                    for st in seg.spiketrains:
                        # Get cell_idx from annotations or parse from name
                        neuron_id = st.annotations.get("cell_idx")
                        if neuron_id is None:
                            # Fallback: parse from "pop_name_cellN_spikes" format
                            if "_cell" in st.name:
                                cell_part = st.name.split("_cell")[-1]
                                # Remove any suffix like "_spikes"
                                neuron_id = int(cell_part.split("_")[0])
                            else:
                                neuron_id = int(st.name)
                        spike_times = st.times.rescale("ms").magnitude
                        times_list.extend(spike_times)
                        ids_list.extend([neuron_id] * len(spike_times))

                    spike_data_arrays[pop_name] = {
                        "times": np.array(times_list),
                        "ids": np.array(ids_list),
                    }

                    if self.verbose:
                        print(
                            f"{pop_name}: {len(times_list)} spikes from {len(seg.spiketrains)} neurons"
                        )
    else:
        # Use manually recorded spike data (legacy)
        for pop_name, data in self.spike_data.items():
            spike_data_arrays[pop_name] = {
                "times": np.array(data["times"]),
                "ids": np.array(data["ids"]),
            }

    joblib.dump(spike_data_arrays, spike_filename, compress=3)

    # Save metadata
    metadata = {
        "total_chunks": self.chunk_id,
        "chunk_duration__ms": self.chunk_duration__ms,
        "recording_config": self.recording_config,
    }
    joblib.dump(metadata, self.save_path / "metadata.pkl")

    if self.verbose:
        print("\nContinuous saving complete:")
        print(f"\tTotal chunks saved: {self.chunk_id}")
        print(f"\tSpike data saved: {spike_filename}")
        print(f"\tPopulations with spikes: {list(spike_data_arrays.keys())}")
        print(f"\tAll data in: {self.save_path}")

convert_chunks_to_neo 

convert_chunks_to_neo(save_path: Path, duration__ms: Optional[float] = None, spike_data_file: Optional[Path] = None, verbose: bool = True) -> Block

Load chunks and convert to NEO Block format (compatible with SimulationRunner output).

This function creates a NEO Block that's identical in structure to what SimulationRunner.run() would return, making it compatible with existing analysis code.

Parameters:

Name Type Description Default
save_path Path

Directory where chunks were saved

 required
duration__ms float

Total simulation duration in ms (if None, inferred from data)

 None
spike_data_file Path

Path to SimulationRunner spike results file (e.g., 'watanabe__spikes_only.pkl') If provided, spike data will be loaded from this NEO Block instead of chunks

 None
verbose bool

If True, display progress bars and status messages. Set to False to disable.

 True

Returns:

Type Description
Block

NEO Block containing spike trains and analog signals
Source code in myogen/utils/continuous_saver.py 
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
def convert_chunks_to_neo(
    save_path: Path, duration__ms: Optional[float] = None, spike_data_file: Optional[Path] = None, verbose: bool = True
) -> Block:
    """
    Load chunks and convert to NEO Block format (compatible with SimulationRunner output).

    This function creates a NEO Block that's identical in structure to what
    SimulationRunner.run() would return, making it compatible with existing
    analysis code.

    Parameters
    ----------
    save_path : Path
        Directory where chunks were saved
    duration__ms : float, optional
        Total simulation duration in ms (if None, inferred from data)
    spike_data_file : Path, optional
        Path to SimulationRunner spike results file (e.g., 'watanabe__spikes_only.pkl')
        If provided, spike data will be loaded from this NEO Block instead of chunks
    verbose : bool, default=True
        If True, display progress bars and status messages. Set to False to disable.

    Returns
    -------
    Block
        NEO Block containing spike trains and analog signals
    """
    save_path = Path(save_path)

    if verbose:
        print("Converting chunks to NEO Block format...")

    # Load metadata
    metadata = joblib.load(save_path / "metadata.pkl")
    total_chunks = metadata["total_chunks"]
    timestep__ms = None

    # Load spike data - either from external file or from chunks
    if spike_data_file is not None:
        # Load spike data from SimulationRunner results (NEO Block)
        if verbose:
            print(f"\tLoading spike data from: {spike_data_file}")
        spike_results = joblib.load(spike_data_file)
        use_neo_spikes = True
    else:
        # Load spike data from chunks (legacy format)
        spike_filename = save_path / "spikes.pkl"
        if spike_filename.exists():
            if verbose:
                print(f"  Loading spike data from: {spike_filename}")
            spike_data = joblib.load(spike_filename)
            use_neo_spikes = False
        else:
            if verbose:
                print("\tWarning: No spike data found")
            spike_data = {}
            use_neo_spikes = False

    # Load first chunk to get timestep info
    first_chunk = joblib.load(save_path / "chunk_0000.pkl")
    timestep__ms = first_chunk["timestep__ms"]

    # Infer duration if not provided
    if duration__ms is None:
        last_chunk = joblib.load(save_path / f"chunk_{total_chunks - 1:04d}.pkl")
        duration__ms = last_chunk["time_end"]

    if verbose:
        print(f"\tDuration: {duration__ms} ms")
        print(f"\tTimestep: {timestep__ms} ms")
        print(f"\tTotal chunks: {total_chunks}")

    # Create NEO Block
    block = Block()

    # Add spike data for each population
    if verbose:
        print("\tAdding spike trains...")

    if use_neo_spikes:
        # Use spike data from SimulationRunner NEO Block
        for seg in spike_results.segments:
            if len(seg.spiketrains) > 0:
                # Create new segment with spike trains
                new_segment = Segment(name=seg.name)
                for st in seg.spiketrains:
                    new_segment.spiketrains.append(st)
                block.segments.append(new_segment)
                if verbose:
                    print(f"\t{seg.name}: {len(new_segment.spiketrains)} spike trains")
    else:
        # Use spike data from chunks (legacy format)
        for pop_name, spikes in spike_data.items():
            segment = Segment(name=pop_name)

            spike_times = spikes["times"]
            spike_ids = spikes["ids"]

            # Create spike trains for each neuron
            unique_ids = sorted(np.unique(spike_ids))
            for spike_id in tqdm(
                unique_ids, desc=f"\tCreating {pop_name} spike trains", leave=False, disable=not verbose
            ):
                times_for_id = spike_times[spike_ids == spike_id]

                # Filter out spike times that exceed duration due to floating-point precision
                # Keep only spikes strictly less than t_stop
                times_for_id = times_for_id[times_for_id < duration__ms]

                if len(times_for_id) > 0:
                    segment.spiketrains.append(
                        SpikeTrain(
                            name=str(int(spike_id)),
                            times=(times_for_id * pq.ms).rescale(pq.s),
                            t_start=0.0 * pq.s,
                            t_stop=(duration__ms * pq.ms).rescale(pq.s),
                            sampling_rate=(1.0 / (timestep__ms * pq.ms)).rescale(pq.Hz),
                        )
                    )

            block.segments.append(segment)
            if verbose:
                print(f"\t{pop_name}: {len(segment.spiketrains)} spike trains")

    # Add membrane potential data by loading and combining chunks
    if verbose:
        print(f"\tLoading and combining membrane data from {total_chunks} chunks...")

    # Determine which populations have membrane recordings
    first_chunk_membrane = first_chunk["membrane_data"]

    for pop_name in first_chunk_membrane.keys():
        # Find or create segment for this population
        segment = None
        for seg in block.segments:
            if seg.name == pop_name:
                segment = seg
                break

        if segment is None:
            segment = Segment(name=pop_name)
            block.segments.append(segment)

        # Get all cell indices for this population
        cell_indices = sorted(first_chunk_membrane[pop_name].keys())

        if verbose:
            print(f"\t{pop_name}: Combining {len(cell_indices)} neurons from {total_chunks} chunks...")

        # OPTIMIZED: Load each chunk once and extract all neurons
        # This reduces file reads from (neurons × chunks) to just (chunks)
        # For 400 neurons × 36 chunks: 14,400 reads → 36 reads (400x faster!)

        # Initialize storage for each neuron
        neuron_data = {cell_idx: [] for cell_idx in cell_indices}

        # Load chunks once and distribute data to neurons
        for chunk_id in tqdm(
            range(total_chunks), desc=f"\tLoading {pop_name} chunks", unit="chunk", disable=not verbose
        ):
            chunk = joblib.load(save_path / f"chunk_{chunk_id:04d}.pkl")

            if pop_name in chunk["membrane_data"]:
                for cell_idx in cell_indices:
                    if cell_idx in chunk["membrane_data"][pop_name]:
                        neuron_data[cell_idx].append(chunk["membrane_data"][pop_name][cell_idx])

        # Concatenate data for each neuron and create AnalogSignals
        if verbose:
            print(f"\t{pop_name}: Creating analog signals...")
        for cell_idx in tqdm(
            cell_indices, desc=f"\tCreating {pop_name} signals", leave=False, unit="signal", disable=not verbose
        ):
            if neuron_data[cell_idx]:
                combined_voltage = np.concatenate(neuron_data[cell_idx])

                segment.analogsignals.append(
                    AnalogSignal(
                        name=str(cell_idx),
                        sampling_period=(timestep__ms * pq.ms).rescale(pq.s),
                        signal=combined_voltage * pq.mV,
                    )
                )

        if verbose:
            print(f"\t{pop_name}: {len(segment.analogsignals)} analog signals created")

    # Add metadata annotations
    block.annotations["time__ms"] = duration__ms
    block.annotations["timestep__ms"] = timestep__ms
    block.annotations["temperature__celsius"] = 36.0  # Default from SimulationRunner

    if verbose:
        print("\nNEO Block created successfully")
        print(f"\tTotal segments: {len(block.segments)}")

    return block

NWB export

Note

NWB export requires optional dependencies: pip install myogen[nwb].

export_to_nwb 

export_to_nwb(block: Block, filepath: str | Path, session_description: str = 'MyoGen simulation', identifier: str | None = None, session_start_time: datetime | None = None, experimenter: str | list[str] | None = None, institution: str | None = None, lab: str | None = None, experiment_description: str | None = None, keywords: list[str] | None = None, subject_id: str | None = None, species: str = 'Homo sapiens', age: str | None = None, sex: str | None = None, subject_description: str | None = None, **kwargs) -> Path

Export a Neo Block to NWB format.

This function uses Neo's NWBIO to write simulation data to an NWB file. The Block should contain AnalogSignals with grid annotations (created via create_grid_signal) for electrode array data.

Parameters:

Name Type Description Default
block Block

Neo Block containing simulation data. Can be spike trains, EMG, or MUAP data from MyoGen simulations.

 required
filepath str or Path

Output file path. Should end with '.nwb'.

 required
session_description str

Description of the simulation session.

 "MyoGen simulation"
identifier str

Unique identifier for this NWB file. If None, a UUID is generated.

 None
session_start_time datetime

Start time of the session. If None, current time is used.

 None
experimenter str or list[str]

Name(s) of experimenter(s).

 None
institution str

Institution where the simulation was performed.

 None
lab str

Lab where the simulation was performed.

 None
experiment_description str

Description of the experiment/simulation.

 None
keywords list[str]

Keywords describing the data.

 None
subject_id str

Subject identifier (recommended for DANDI).

 None
species str

Species of the subject. Use Latin binomial (e.g., "Homo sapiens", "Mus musculus"). For simulations, defaults to human.

 "Homo sapiens"
age str

Age of subject in ISO 8601 duration format (e.g., "P30Y" for 30 years).

 None
sex str

Sex of subject. One of: "M", "F", "U" (unknown), "O" (other).

 None
subject_description str

Description of the subject.

 None
**kwargs

Additional keyword arguments passed to NWBIO.

 {}

Returns:

Type Description
Path

Path to the created NWB file.

Examples:

>>> from myogen.utils.nwb import export_to_nwb
>>>
>>> # Export spike trains to NWB
>>> export_to_nwb(
...     spike_train__Block,
...     "simulation_spikes.nwb",
...     session_description="Motor neuron pool simulation",
...     institution="My University",
... )
>>>
>>> # Export surface EMG to NWB
>>> export_to_nwb(
...     surface_emg__Block,
...     "simulation_emg.nwb",
...     session_description="Surface EMG simulation",
...     experimenter="John Doe",
... )
Notes

For electrode array data (surface EMG, MUAPs), the grid structure is preserved via electrode_positions in annotations, which map to NWB's electrode table.

See Also

create_grid_signal : Create grid-annotated AnalogSignals validate_nwb : Validate NWB file with NWBInspector

Source code in myogen/utils/nwb.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
@beartowertype
def export_to_nwb(
    block: Block,
    filepath: str | Path,
    session_description: str = "MyoGen simulation",
    identifier: str | None = None,
    session_start_time: datetime | None = None,
    experimenter: str | list[str] | None = None,
    institution: str | None = None,
    lab: str | None = None,
    experiment_description: str | None = None,
    keywords: list[str] | None = None,
    # Subject metadata (important for DANDI compliance)
    subject_id: str | None = None,
    species: str = "Homo sapiens",
    age: str | None = None,
    sex: str | None = None,
    subject_description: str | None = None,
    **kwargs,
) -> Path:
    """
    Export a Neo Block to NWB format.

    This function uses Neo's NWBIO to write simulation data to an NWB file.
    The Block should contain AnalogSignals with grid annotations (created
    via create_grid_signal) for electrode array data.

    Parameters
    ----------
    block : Block
        Neo Block containing simulation data. Can be spike trains, EMG,
        or MUAP data from MyoGen simulations.
    filepath : str or Path
        Output file path. Should end with '.nwb'.
    session_description : str, default="MyoGen simulation"
        Description of the simulation session.
    identifier : str, optional
        Unique identifier for this NWB file. If None, a UUID is generated.
    session_start_time : datetime, optional
        Start time of the session. If None, current time is used.
    experimenter : str or list[str], optional
        Name(s) of experimenter(s).
    institution : str, optional
        Institution where the simulation was performed.
    lab : str, optional
        Lab where the simulation was performed.
    experiment_description : str, optional
        Description of the experiment/simulation.
    keywords : list[str], optional
        Keywords describing the data.
    subject_id : str, optional
        Subject identifier (recommended for DANDI).
    species : str, default="Homo sapiens"
        Species of the subject. Use Latin binomial (e.g., "Homo sapiens",
        "Mus musculus"). For simulations, defaults to human.
    age : str, optional
        Age of subject in ISO 8601 duration format (e.g., "P30Y" for 30 years).
    sex : str, optional
        Sex of subject. One of: "M", "F", "U" (unknown), "O" (other).
    subject_description : str, optional
        Description of the subject.
    **kwargs
        Additional keyword arguments passed to NWBIO.

    Returns
    -------
    Path
        Path to the created NWB file.

    Examples
    --------
    >>> from myogen.utils.nwb import export_to_nwb
    >>>
    >>> # Export spike trains to NWB
    >>> export_to_nwb(
    ...     spike_train__Block,
    ...     "simulation_spikes.nwb",
    ...     session_description="Motor neuron pool simulation",
    ...     institution="My University",
    ... )
    >>>
    >>> # Export surface EMG to NWB
    >>> export_to_nwb(
    ...     surface_emg__Block,
    ...     "simulation_emg.nwb",
    ...     session_description="Surface EMG simulation",
    ...     experimenter="John Doe",
    ... )

    Notes
    -----
    For electrode array data (surface EMG, MUAPs), the grid structure is
    preserved via electrode_positions in annotations, which map to NWB's
    electrode table.

    See Also
    --------
    create_grid_signal : Create grid-annotated AnalogSignals
    validate_nwb : Validate NWB file with NWBInspector
    """
    _check_nwb_available()

    filepath = Path(filepath)
    if not filepath.suffix == ".nwb":
        filepath = filepath.with_suffix(".nwb")

    # Generate defaults
    if identifier is None:
        # Derive a stable identifier from the current MyoGen seed so the same
        # simulation produces the same NWB identifier (DANDI-friendly). Falls
        # back to a fixed namespace UUID if the seed module is unavailable.
        try:
            from myogen import get_random_seed

            identifier = str(
                uuid.uuid5(uuid.NAMESPACE_OID, f"myogen-{get_random_seed()}")
            )
        except Exception:
            identifier = str(uuid.uuid5(uuid.NAMESPACE_OID, "myogen"))
    if session_start_time is None:
        # Default to a fixed epoch so two runs produce byte-identical NWB
        # files; warn so callers know to pass a real time for live recordings.
        _warnings.warn(
            "export_to_nwb called without an explicit session_start_time; "
            "defaulting to 1970-01-01T00:00:00Z so the export is reproducible. "
            "Pass session_start_time=datetime.now(tz=...) (or similar) for "
            "real recordings.",
            UserWarning,
            stacklevel=2,
        )
        session_start_time = datetime(1970, 1, 1, tzinfo=timezone.utc)
    if keywords is None:
        keywords = ["MyoGen", "simulation", "EMG", "motor unit"]

    # Build metadata dict
    nwb_metadata = {
        "session_description": session_description,
        "identifier": identifier,
        "session_start_time": session_start_time,
    }
    if experimenter is not None:
        nwb_metadata["experimenter"] = (
            [experimenter] if isinstance(experimenter, str) else experimenter
        )
    if institution is not None:
        nwb_metadata["institution"] = institution
    if lab is not None:
        nwb_metadata["lab"] = lab
    if experiment_description is not None:
        nwb_metadata["experiment_description"] = experiment_description
    if keywords:
        nwb_metadata["keywords"] = keywords

    # Merge with any additional kwargs
    nwb_metadata.update(kwargs)

    # Write to NWB
    writer = NWBIO(str(filepath), mode="w", **nwb_metadata)
    writer.write(block)

    # Add subject metadata if provided (important for DANDI compliance)
    # Neo's NWBIO doesn't support subject directly, so we add it post-hoc
    if subject_id is not None or subject_description is not None:
        from pynwb import NWBHDF5IO
        from pynwb.file import Subject

        with NWBHDF5IO(str(filepath), mode="r+") as io:
            nwbfile = io.read()
            nwbfile.subject = Subject(
                subject_id=subject_id or "simulation",
                species=species,
                age=age,
                sex=sex,
                description=subject_description or "Simulated subject",
            )
            io.write(nwbfile)

    return filepath

export_simulation_to_nwb 

export_simulation_to_nwb(filepath: str | Path, spike_train__Block: Block | None = None, surface_emg__Block: Block | None = None, surface_muap__Block: Block | None = None, intramuscular_emg__Block: Block | None = None, intramuscular_muap__Block: Block | None = None, session_description: str = 'MyoGen neuromuscular simulation', identifier: str | None = None, session_start_time: datetime | None = None, **kwargs) -> Path

Export all simulation data to a single NWB file.

This is a convenience function that combines multiple Neo Blocks (spike trains, EMG, MUAPs) into a single NWB file.

Parameters:

Name Type Description Default
filepath str or Path

Output file path.

 required
spike_train__Block SPIKE_TRAIN__Block

Block containing spike train data.

 None
surface_emg__Block SURFACE_EMG__Block

Block containing surface EMG data.

 None
surface_muap__Block SURFACE_MUAP__Block

Block containing surface MUAP templates.

 None
intramuscular_emg__Block INTRAMUSCULAR_EMG__Block

Block containing intramuscular EMG data.

 None
intramuscular_muap__Block INTRAMUSCULAR_MUAP__Block

Block containing intramuscular MUAP templates.

 None
session_description str

Description of the simulation session.

 "MyoGen neuromuscular simulation"
identifier str

Unique identifier for this NWB file.

 None
session_start_time datetime

Start time of the session.

 None
**kwargs

Additional metadata passed to export_to_nwb.

 {}

Returns:

Type Description
Path

Path to the created NWB file.

Examples:

>>> from myogen.utils.nwb import export_simulation_to_nwb
>>>
>>> # Export complete simulation
>>> export_simulation_to_nwb(
...     "full_simulation.nwb",
...     spike_train__Block=simulation.get_spike_train__Block(),
...     surface_emg__Block=surface_emg.surface_emg__Block,
...     session_description="Biceps brachii simulation",
...     institution="University",
... )
Source code in myogen/utils/nwb.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
@beartowertype
def export_simulation_to_nwb(
    filepath: str | Path,
    spike_train__Block: Block | None = None,
    surface_emg__Block: Block | None = None,
    surface_muap__Block: Block | None = None,
    intramuscular_emg__Block: Block | None = None,
    intramuscular_muap__Block: Block | None = None,
    session_description: str = "MyoGen neuromuscular simulation",
    identifier: str | None = None,
    session_start_time: datetime | None = None,
    **kwargs,
) -> Path:
    """
    Export all simulation data to a single NWB file.

    This is a convenience function that combines multiple Neo Blocks
    (spike trains, EMG, MUAPs) into a single NWB file.

    Parameters
    ----------
    filepath : str or Path
        Output file path.
    spike_train__Block : SPIKE_TRAIN__Block, optional
        Block containing spike train data.
    surface_emg__Block : SURFACE_EMG__Block, optional
        Block containing surface EMG data.
    surface_muap__Block : SURFACE_MUAP__Block, optional
        Block containing surface MUAP templates.
    intramuscular_emg__Block : INTRAMUSCULAR_EMG__Block, optional
        Block containing intramuscular EMG data.
    intramuscular_muap__Block : INTRAMUSCULAR_MUAP__Block, optional
        Block containing intramuscular MUAP templates.
    session_description : str, default="MyoGen neuromuscular simulation"
        Description of the simulation session.
    identifier : str, optional
        Unique identifier for this NWB file.
    session_start_time : datetime, optional
        Start time of the session.
    **kwargs
        Additional metadata passed to export_to_nwb.

    Returns
    -------
    Path
        Path to the created NWB file.

    Examples
    --------
    >>> from myogen.utils.nwb import export_simulation_to_nwb
    >>>
    >>> # Export complete simulation
    >>> export_simulation_to_nwb(
    ...     "full_simulation.nwb",
    ...     spike_train__Block=simulation.get_spike_train__Block(),
    ...     surface_emg__Block=surface_emg.surface_emg__Block,
    ...     session_description="Biceps brachii simulation",
    ...     institution="University",
    ... )
    """
    _check_nwb_available()

    # Combine all blocks into one
    combined_block = Block(name="MyoGen_Simulation")

    # Add annotations to identify data types
    if spike_train__Block is not None:
        for segment in spike_train__Block.segments:
            segment.annotate(myogen_data_type="spike_train")
            combined_block.segments.append(segment)

    if surface_emg__Block is not None:
        for group in surface_emg__Block.groups:
            group.annotate(myogen_data_type="surface_emg")
            combined_block.groups.append(group)

    if surface_muap__Block is not None:
        for group in surface_muap__Block.groups:
            group.annotate(myogen_data_type="surface_muap")
            combined_block.groups.append(group)

    if intramuscular_emg__Block is not None:
        for segment in intramuscular_emg__Block.segments:
            segment.annotate(myogen_data_type="intramuscular_emg")
            combined_block.segments.append(segment)

    if intramuscular_muap__Block is not None:
        for segment in intramuscular_muap__Block.segments:
            segment.annotate(myogen_data_type="intramuscular_muap")
            combined_block.segments.append(segment)

    return export_to_nwb(
        combined_block,
        filepath,
        session_description=session_description,
        identifier=identifier,
        session_start_time=session_start_time,
        **kwargs,
    )

validate_nwb 

validate_nwb(filepath: str | Path, verbose: bool = True) -> bool

Validate an NWB file using NWBInspector.

Parameters:

Name Type Description Default
filepath str or Path

Path to the NWB file to validate.

 required
verbose bool

If True, print validation results.

 True

Returns:

Type Description
bool

True if validation passed with no errors, False otherwise.

Examples:

>>> from myogen.utils.nwb import validate_nwb
>>>
>>> is_valid = validate_nwb("simulation.nwb")
>>> if is_valid:
...     print("File is valid!")
Source code in myogen/utils/nwb.py 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
@beartowertype
def validate_nwb(filepath: str | Path, verbose: bool = True) -> bool:
    """
    Validate an NWB file using NWBInspector.

    Parameters
    ----------
    filepath : str or Path
        Path to the NWB file to validate.
    verbose : bool, default=True
        If True, print validation results.

    Returns
    -------
    bool
        True if validation passed with no errors, False otherwise.

    Examples
    --------
    >>> from myogen.utils.nwb import validate_nwb
    >>>
    >>> is_valid = validate_nwb("simulation.nwb")
    >>> if is_valid:
    ...     print("File is valid!")
    """
    try:
        from nwbinspector import inspect_nwbfile
        from nwbinspector.inspector_tools import format_messages
    except ImportError:
        if verbose:
            print(
                "NWBInspector not installed. Install with: pip install nwbinspector"
            )
        return True  # Can't validate without inspector

    filepath = Path(filepath)
    if not filepath.exists():
        raise FileNotFoundError(f"NWB file not found: {filepath}")

    # Run inspection
    messages = list(inspect_nwbfile(nwbfile_path=str(filepath)))

    # Filter by severity
    errors = [m for m in messages if m.importance.name == "CRITICAL"]
    warnings = [m for m in messages if m.importance.name in ("ERROR", "WARNING")]

    if verbose:
        if not messages:
            print(f"✓ {filepath.name}: No issues found")
        else:
            print(f"\n{filepath.name} validation results:")
            print(format_messages(messages))

    return len(errors) == 0