simulation

Simulation execution and results handling.

toric_spines_sim.simulation

Simulation execution and results handling.

SimulationResults(results_dict)

Container for simulation results with analysis methods.

This class wraps the raw simulation results and provides convenient methods for analyzing voltage traces, including integration by segment tag with optional surface-area weighting.

Attributes: voltage_traces: TsdFrame with voltage data (rows=time, columns=probes) input_events: TsGroup mapping stream indices to event timestamps synapses: Dict of synapse specifications gap_junctions: List of gap junction specifications record_points: Dict mapping probe labels to (x, y, z) coordinates cell: arbor.cable_cell object morphology: arbor.morphology object segment_tree: arbor.segment_tree object decor: arbor.decor object labels: arbor.label_dict object cvp: arbor.cv_policy object swc_filepath: Path to SWC morphology file synpts_filepath: Path to synapse points file

Initialize from a simulation results dictionary.

Source code in toric_spines_sim/simulation/results.py

def __init__(self, results_dict: Dict):
    """Initialize from a simulation results dictionary."""
    self.voltage_traces = results_dict["voltage_traces"]
    self.input_events = results_dict["input_events"]
    self.synapses = results_dict["synapses"]
    self.gap_junctions = results_dict["gap_junctions"]
    self.record_points = results_dict["record_points"]
    self.cell = results_dict["cell"]
    self.morphology = results_dict["morphology"]
    self.segment_tree = results_dict["segment_tree"]
    self.decor = results_dict["decor"]
    self.labels = results_dict["labels"]
    self.cvp = results_dict["cvp"]
    self.swc_filepath = results_dict.get("swc_filepath")
    self.synpts_filepath = results_dict.get("synpts_filepath")

    # Cache for segment properties
    self._segment_tags_cache = None
    self._segment_areas_cache = None
    self._probe_to_segment_cache = None

integrate_voltages_by_tag(tags, method='average')

Integrate voltage traces from all segments with specified tag(s).

Args: tags: Single tag or collection of tags to integrate over. method: Integration method: - "average": Simple arithmetic mean of all voltage traces - "surface_weighted": Weighted average by segment surface area

Returns: Pynapple Tsd with time in seconds and integrated voltage in mV.

Example: >>> results = SimulationResults(raw_results) >>> # Get average voltage across all spine segments (tag 3) >>> v_spine = results.integrate_voltages_by_tag(3, method="average") >>> >>> # Get surface-weighted voltage across sink segments (tag 5) >>> v_sink = results.integrate_voltages_by_tag(5, method="surface_weighted") >>> >>> # Combine multiple tags >>> v_combined = results.integrate_voltages_by_tag([3, 5])

Source code in toric_spines_sim/simulation/results.py

def integrate_voltages_by_tag(
    self,
    tags: Union[int, List[int], Set[int]],
    method: Literal["average", "surface_weighted"] = "average",
) -> nap.Tsd:
    """
    Integrate voltage traces from all segments with specified tag(s).

    Args:
        tags: Single tag or collection of tags to integrate over.
        method: Integration method:
            - "average": Simple arithmetic mean of all voltage traces
            - "surface_weighted": Weighted average by segment surface area

    Returns:
        Pynapple Tsd with time in seconds and integrated voltage in mV.

    Example:
        >>> results = SimulationResults(raw_results)
        >>> # Get average voltage across all spine segments (tag 3)
        >>> v_spine = results.integrate_voltages_by_tag(3, method="average")
        >>>
        >>> # Get surface-weighted voltage across sink segments (tag 5)
        >>> v_sink = results.integrate_voltages_by_tag(5, method="surface_weighted")
        >>>
        >>> # Combine multiple tags
        >>> v_combined = results.integrate_voltages_by_tag([3, 5])
    """
    # Normalize tags to a set
    if isinstance(tags, int):
        tags_set = {tags}
    else:
        tags_set = set(tags)

    # Get mappings
    segment_tags = self._get_segment_tags()
    probe_to_segment = self._map_probes_to_segments()

    # Find probes that correspond to segments with the target tags
    matching_probes = []
    for probe_label, seg_idx in probe_to_segment.items():
        if segment_tags.get(seg_idx) in tags_set:
            matching_probes.append((probe_label, seg_idx))

    if not matching_probes:
        raise ValueError(f"No segments found with tag(s) {tags_set}")

    logger.info(
        f"Integrating {len(matching_probes)} voltage traces for tag(s) {tags_set}"
    )

    # Extract matching column names from TsdFrame
    probe_labels = [p[0] for p in matching_probes]

    # Check which probes exist in the voltage_traces
    available_probes = [p for p in probe_labels if p in self.voltage_traces.columns]
    if not available_probes:
        raise ValueError(f"No voltage data found for tag(s) {tags_set}")

    # Get subset of TsdFrame
    voltage_subset = self.voltage_traces[available_probes]

    # Get time array from TsdFrame index
    time = voltage_subset.index.values

    # Integrate based on method
    if method == "average":
        # Simple average using TsdFrame.mean(axis=1)
        integrated_voltage = voltage_subset.mean(axis=1).values
        logger.debug(f"Computed simple average of {len(available_probes)} traces")

    elif method == "surface_weighted":
        # Weighted average by surface area
        segment_areas = self._get_segment_surface_areas()

        # Map probe labels to their segment indices
        probe_to_seg = {p: s for p, s in matching_probes}

        weights = []
        for probe_label in available_probes:
            seg_idx = probe_to_seg[probe_label]
            weights.append(segment_areas[seg_idx])

        weights = np.array(weights)
        weights = weights / weights.sum()  # Normalize

        # Weighted sum: multiply each column by its weight and sum
        weighted_data = voltage_subset.values * weights[np.newaxis, :]
        integrated_voltage = weighted_data.sum(axis=1)

        logger.debug(
            f"Computed surface-weighted average of {len(available_probes)} traces"
        )
    else:
        raise ValueError(
            f"Unknown method '{method}'. Use 'average' or 'surface_weighted'"
        )

    # Return as pynapple Tsd (time in seconds, voltage in mV)
    return nap.Tsd(t=time, d=integrated_voltage, time_units="s")

get_tags()

Get all unique tags present in the morphology.

Source code in toric_spines_sim/simulation/results.py

def get_tags(self) -> Set[int]:
    """Get all unique tags present in the morphology."""
    return set(self._get_segment_tags().values())

get_segments_by_tag(tag)

Get list of segment indices with the specified tag.

Source code in toric_spines_sim/simulation/results.py

def get_segments_by_tag(self, tag: int) -> List[int]:
    """Get list of segment indices with the specified tag."""
    segment_tags = self._get_segment_tags()
    return [seg_idx for seg_idx, seg_tag in segment_tags.items() if seg_tag == tag]

to_dict()

Convert back to a dictionary (for saving).

Source code in toric_spines_sim/simulation/results.py

def to_dict(self) -> Dict:
    """Convert back to a dictionary (for saving)."""
    return {
        "voltage_traces": self.voltage_traces,
        "input_events": self.input_events,
        "synapses": self.synapses,
        "gap_junctions": self.gap_junctions,
        "record_points": self.record_points,
        "cell": self.cell,
        "morphology": self.morphology,
        "segment_tree": self.segment_tree,
        "decor": self.decor,
        "labels": self.labels,
        "cvp": self.cvp,
    }

to_serializable_dict()

Convert to a dictionary containing only serializable data.

Excludes Arbor C++ objects (cell, morphology, segment_tree, decor, labels, cvp) and SynapsePoint/GapJunctionPoint objects which cannot be pickled. These can be reconstructed from the SWC and synapse points files if needed.

Source code in toric_spines_sim/simulation/results.py

def to_serializable_dict(self) -> Dict:
    """Convert to a dictionary containing only serializable data.

    Excludes Arbor C++ objects (cell, morphology, segment_tree, decor, labels, cvp)
    and SynapsePoint/GapJunctionPoint objects which cannot be pickled.
    These can be reconstructed from the SWC and synapse points files if needed.
    """
    # Convert TsdFrame to serializable format
    # Pynapple stores times internally in seconds; retrieve as ms
    times_ms = self.voltage_traces.as_units("ms").index.values
    serializable_voltages = {
        "t": times_ms.tolist(),
        "d": self.voltage_traces.values.tolist(),
        "columns": list(self.voltage_traces.columns),
        "time_units": "ms",
    }

    # Convert TsGroup to serializable format
    serializable_events = {}
    for idx, ts in self.input_events.items():
        label = (
            self.input_events.get_info("label")[idx]
            if "label" in self.input_events.metadata
            else str(idx)
        )
        # Retrieve event times in ms
        times_ms = ts.as_units("ms").index.values
        serializable_events[idx] = {
            "t": times_ms.tolist(),
            "label": label,
            "time_units": "ms",
        }

    return {
        "voltage_traces": serializable_voltages,
        "input_events": serializable_events,
        "record_points": self.record_points,
        "swc_filepath": str(self.swc_filepath) if self.swc_filepath else None,
        "synpts_filepath": (
            str(self.synpts_filepath) if self.synpts_filepath else None
        ),
    }

save(filepath)

Save simulation results to a pickle file.

Note: Only serializable data is saved (voltages, events, record_points, metadata). Arbor C++ objects (cell, morphology, segment_tree, etc.) and SynapsePoint/GapJunctionPoint objects cannot be pickled and are excluded. To access these objects later, you'll need to keep the original SimulationResults instance or rebuild them from the SWC and synapse files.

Args: filepath: Path where to save the results (.pkl file).

Example: >>> results.save("simulations/ts1/results.pkl")

Source code in toric_spines_sim/simulation/results.py

def save(self, filepath: Union[str, Path]) -> None:
    """
    Save simulation results to a pickle file.

    Note: Only serializable data is saved (voltages, events, record_points, metadata).
    Arbor C++ objects (cell, morphology, segment_tree, etc.) and SynapsePoint/GapJunctionPoint
    objects cannot be pickled and are excluded. To access these objects later, you'll need to
    keep the original SimulationResults instance or rebuild them from the SWC and synapse files.

    Args:
        filepath: Path where to save the results (.pkl file).

    Example:
        >>> results.save("simulations/ts1/results.pkl")
    """
    filepath = Path(filepath)
    logger.info(f"Saving simulation results to {filepath}")

    try:
        # Ensure parent directory exists
        filepath.parent.mkdir(parents=True, exist_ok=True)

        # Save only serializable data (Arbor C++ objects cannot be pickled)
        serializable_data = self.to_serializable_dict()

        with open(filepath, "wb") as f:
            dill.dump(serializable_data, f)

        file_size_mb = filepath.stat().st_size / (1024 * 1024)
        logger.info(f"Results saved successfully ({file_size_mb:.2f} MB)")
        logger.warning(
            "Note: Arbor objects and synapse/gap junction objects were not saved as they cannot be serialized. "
            "Use swc_filepath and synpts_filepath to reconstruct if needed."
        )
    except Exception as e:
        logger.error(f"Failed to save results: {e}")
        raise

load(filepath) classmethod

Load simulation results from a pickle file.

Note: Loaded results will not contain Arbor objects (cell, morphology, etc.) or synapse/gap junction objects as these cannot be serialized. Only voltage_traces (TsdFrame), input_events (TsGroup), record_points, and metadata are loaded. Methods that require Arbor objects (like integrate_voltages_by_tag) will not work on loaded results.

Args: filepath: Path to the saved results file (.pkl).

Returns: SimulationResults instance with limited functionality.

Raises: FileNotFoundError: If the file does not exist.

Example: >>> results = SimulationResults.load("simulations/ts1/results.pkl") >>> voltage_trace = results.voltage_traces["probe_seg_0"]

Source code in toric_spines_sim/simulation/results.py

@classmethod
def load(cls, filepath: Union[str, Path]) -> "SimulationResults":
    """
    Load simulation results from a pickle file.

    Note: Loaded results will not contain Arbor objects (cell, morphology, etc.)
    or synapse/gap junction objects as these cannot be serialized. Only voltage_traces
    (TsdFrame), input_events (TsGroup), record_points, and metadata are loaded.
    Methods that require Arbor objects (like integrate_voltages_by_tag) will not work
    on loaded results.

    Args:
        filepath: Path to the saved results file (.pkl).

    Returns:
        SimulationResults instance with limited functionality.

    Raises:
        FileNotFoundError: If the file does not exist.

    Example:
        >>> results = SimulationResults.load("simulations/ts1/results.pkl")
        >>> voltage_trace = results.voltage_traces["probe_seg_0"]
    """
    filepath = Path(filepath)

    if not filepath.exists():
        raise FileNotFoundError(f"Simulation results file not found: {filepath}")

    logger.info(f"Loading simulation results from {filepath}")

    try:
        with open(filepath, "rb") as f:
            results_dict = dill.load(f)

        file_size_mb = filepath.stat().st_size / (1024 * 1024)
        logger.info(f"Results loaded successfully ({file_size_mb:.2f} MB)")

        # Reconstruct TsdFrame from serializable format
        voltage_data = results_dict.get("voltage_traces", {})
        if isinstance(voltage_data, dict) and "t" in voltage_data:
            # Reconstruct from serialized format
            voltage_traces = nap.TsdFrame(
                t=voltage_data["t"],
                d=voltage_data["d"],
                time_units=voltage_data.get("time_units", "ms"),
                columns=voltage_data["columns"],
            )
        else:
            # Legacy format - empty TsdFrame
            voltage_traces = nap.TsdFrame(t=[], d=[], time_units="ms")

        # Reconstruct TsGroup from serializable format
        events_data = results_dict.get("input_events", {})
        if isinstance(events_data, dict):
            ts_dict = {}
            labels = []
            for idx, event_info in events_data.items():
                if isinstance(event_info, dict):
                    ts_dict[int(idx)] = nap.Ts(
                        t=event_info["t"],
                        time_units=event_info.get("time_units", "ms"),
                    )
                    labels.append(event_info.get("label", str(idx)))
                else:
                    # Legacy format: just a list of times
                    ts_dict[int(idx)] = nap.Ts(t=event_info, time_units="ms")
                    labels.append(str(idx))
            input_events = nap.TsGroup(ts_dict, label=labels)
        else:
            # Legacy format - empty TsGroup
            input_events = nap.TsGroup({})

        # Add None placeholders for Arbor objects that weren't saved
        full_dict = {
            "voltage_traces": voltage_traces,
            "input_events": input_events,
            "synapses": {},
            "gap_junctions": {},
            "record_points": results_dict.get("record_points", {}),
            "cell": None,
            "morphology": None,
            "segment_tree": None,
            "decor": None,
            "labels": None,
            "cvp": None,
            "swc_filepath": results_dict.get("swc_filepath"),
            "synpts_filepath": results_dict.get("synpts_filepath"),
        }

        # Log some basic info about the loaded results
        logger.info(f"Loaded voltage TsdFrame with shape {voltage_traces.shape}")
        logger.info(f"Loaded events TsGroup with {len(input_events)} streams")
        if "swc_filepath" in results_dict:
            logger.info(f"SWC file: {results_dict['swc_filepath']}")
        if "synpts_filepath" in results_dict:
            logger.info(f"Synapse points file: {results_dict['synpts_filepath']}")
        logger.warning(
            "Note: Arbor objects and synapse/gap junction objects were not loaded (not available in saved file)"
        )

        return cls(full_dict)

    except Exception as e:
        logger.error(f"Failed to load results: {e}")
        raise

TSSimulator(swc_filepath, synpts_filepath, events, parameters, record_points='all')

Run an Arbor simulation of a toric-spine SWC with AMPA synapses.

All configuration is provided at construction time. Intermediate pipeline objects (synapses, cell, recipe, etc.) are computed on first use, then reused for the lifetime of the instance. run() may be called multiple times; each call reuses cached pipeline objects and creates a fresh Arbor simulation.

Synapses are always built as AMPA from the points file (syn_0, syn_1, … in file order). For other receptor types, build a TSModel / SynapsePopulation yourself.

Event channels are mapped by TsGroup index to that synapse order. Generators emit axon-order channels (A0S0, …). If axon assignments are not already point-file order, call remap_axon_channel_events_to_synapses before passing events here.

When using stochastic event generators, pass parameters["seed"] to the generator so that repeated run() calls with the same instance are deterministic.

Examples:

>>> from toric_spines_sim.simulation import TSSimulator
>>> from toric_spines_sim.simulation import make_default_parameter_bank
>>> from toric_spines_sim.events import StochasticEventGenerator, FlatRateCurve
>>>
>>> parameter_bank = make_default_parameter_bank()
>>> parameters = parameter_bank.sample()
>>> events = StochasticEventGenerator(
...     rate_curves=[FlatRateCurve(rate_hz=50.0)],
...     n_synapses_per_axon=[25],
...     T_ms=parameters["T_ms"],
...     seed=int(parameters["seed"]),
... ).generate()
>>> sim = TSSimulator(
...     "data/swc/microns/TS1_wsink_r10um.swc",
...     "data/pointsets/microns/TS1_synpts.txt",
...     events,
...     parameters,
... )
>>> results = sim.run()

Source code in toric_spines_sim/simulation/core.py

def __init__(
    self,
    swc_filepath: Union[str, Path],
    synpts_filepath: Union[str, Path],
    events: nap.TsGroup,
    parameters: ParameterSet,
    record_points: Union[
        Literal["all"], Dict[str, Tuple[float, float, float]]
    ] = "all",
):
    self._swc_filepath = Path(swc_filepath)
    self._synpts_filepath = Path(synpts_filepath)
    self._events = events
    self._parameters = parameters
    self._record_points_spec = record_points

    # Pipeline caches (computed on first use, then reused)
    self._synapses: Optional[Dict] = None
    self._gap_junctions: Optional[Dict] = None
    self._record_points_resolved: Optional[
        Dict[str, Tuple[float, float, float]]
    ] = None
    self._build_cell_results: Optional[Dict] = None
    self._recipe: Optional[TSRecipe] = None

    logger.info(
        "Initialized TSSimulator: swc=%s, synpts=%s, events=%d streams, record_points=%s",
        self._swc_filepath,
        self._synpts_filepath,
        len(events),
        record_points if record_points != "all" else "all",
    )

swc_filepath property

Path to the SWC morphology file.

synpts_filepath property

Path to the synapse points file.

parameters property

Sampled simulation parameters.

record_points_spec property

Recording point specification passed at construction.

synapses property

Synapse specifications (builds if needed).

gap_junctions property

Gap junction specifications (builds synapses if needed).

record_points property

Resolved record points (builds if needed).

events property

Input event timestamps.

cell property

The built Arbor cable_cell (builds if needed).

morphology property

The cell morphology (builds cell if needed).

segment_tree property

The segment tree (builds cell if needed).

decor property

The cell decor (builds cell if needed).

labels property

The label dictionary (builds cell if needed).

cvp property

The control volume policy (builds cell if needed).

recipe property

The Arbor recipe (builds if needed).

build_synapses()

Build synapses and gap junctions from morphology files.

Returns: Dictionary of synapse specifications.

Source code in toric_spines_sim/simulation/core.py

def build_synapses(self) -> Dict:
    """Build synapses and gap junctions from morphology files.

    Returns:
        Dictionary of synapse specifications.
    """
    if self._synapses is not None:
        return self._synapses

    logger.info("Building synapses and gap junctions...")
    self._synapses = SynapsePopulation.from_file(
        self._synpts_filepath,
        model="ampa",
        global_parameters=self._parameters,
    ).synapses
    self._gap_junctions = prepare_gap_junctions(
        self._swc_filepath, parameters=self._parameters
    )
    logger.info(
        "Built %d synapses and %d gap junctions",
        len(self._synapses),
        len(self._gap_junctions),
    )

    return self._synapses

build_record_points()

Resolve record points (either "all" or explicit dict).

Returns: Dictionary mapping probe labels to (x, y, z) coordinates.

Source code in toric_spines_sim/simulation/core.py

def build_record_points(self) -> Dict[str, Tuple[float, float, float]]:
    """Resolve record points (either "all" or explicit dict).

    Returns:
        Dictionary mapping probe labels to (x, y, z) coordinates.
    """
    if self._record_points_resolved is not None:
        return self._record_points_resolved

    if self._record_points_spec == "all":
        logger.info("Resolving record points from all segment centers")
        self._record_points_resolved = get_center_coordinates_for_all_segments(
            self._swc_filepath
        )
        logger.info(
            "Resolved %d recording points",
            len(self._record_points_resolved),
        )
    else:
        self._record_points_resolved = self._record_points_spec
        logger.info(
            "Using %d specified recording points",
            len(self._record_points_resolved),
        )

    return self._record_points_resolved

build_events()

Return input events.

Returns: TsGroup mapping stream indices to Ts objects (timestamps in ms).

Source code in toric_spines_sim/simulation/core.py

def build_events(self) -> nap.TsGroup:
    """Return input events.

    Returns:
        TsGroup mapping stream indices to Ts objects (timestamps in ms).
    """
    return self._events

build_cell()

Build the Arbor cable cell.

Returns: The built Arbor cable_cell.

Source code in toric_spines_sim/simulation/core.py

def build_cell(self) -> A.cable_cell:
    """Build the Arbor cable cell.

    Returns:
        The built Arbor cable_cell.
    """
    if self._build_cell_results is not None:
        return self._build_cell_results["cell"]

    synapses = self.build_synapses()
    gap_junctions = self._gap_junctions
    record_points = self.build_record_points()

    logger.info("Building cell morphology...")
    tsm = TSModel(
        swc_path=self._swc_filepath,
        synapses=synapses,
        gap_junctions=gap_junctions,
        record_points=record_points,
        parameters=self._parameters,
    )
    self._build_cell_results = tsm.build_cell()
    logger.info("Cell built successfully")

    return self._build_cell_results["cell"]

build_recipe()

Build the Arbor recipe.

Returns: The built TSRecipe.

Source code in toric_spines_sim/simulation/core.py

def build_recipe(self) -> TSRecipe:
    """Build the Arbor recipe.

    Returns:
        The built TSRecipe.
    """
    if self._recipe is not None:
        return self._recipe

    cell = self.build_cell()
    synapses = self.build_synapses()
    gap_junctions = self._gap_junctions
    record_points = self.build_record_points()
    events = self.build_events()

    logger.info("Creating Arbor recipe...")
    self._recipe = TSRecipe(
        cell,
        synapses=synapses,
        gap_junctions=gap_junctions,
        record_points=record_points,
        events=events,
        parameters=self._parameters,
        custom_catalogue=self._build_cell_results["custom_catalogue"],
    )

    return self._recipe

write_cell(filename)

Write the cable cell to file using Arbor's write_component.

Args: filename: Path to the output file (should have .acc extension)

Source code in toric_spines_sim/simulation/core.py

def write_cell(self, filename: str) -> None:
    """Write the cable cell to file using Arbor's write_component.

    Args:
        filename: Path to the output file (should have .acc extension)
    """
    A.write_component(self.cell, filename)
    logger.info("Wrote cable cell to %s", filename)

run()

Run the simulation and return results.

May be called multiple times on the same instance. Pipeline objects are reused; a fresh Arbor simulation is created on each call.

Returns: SimulationResults containing all simulation data.

Source code in toric_spines_sim/simulation/core.py

def run(self) -> SimulationResults:
    """Run the simulation and return results.

    May be called multiple times on the same instance. Pipeline objects
    are reused; a fresh Arbor simulation is created on each call.

    Returns:
        SimulationResults containing all simulation data.
    """
    recipe = self.build_recipe()
    record_points = self.build_record_points()
    events = self.build_events()

    logger.info("Setting up Arbor simulation...")
    ctx = A.context()
    dec = A.partition_load_balance(recipe, ctx)
    sim = A.simulation(recipe, ctx, dec)
    logger.info("Simulation context: %s", ctx)

    logger.info("Setting up %d voltage probes...", len(record_points))
    probe_handles = {}
    dt_record_ms = self._parameters["dt_record_ms"]
    for probe_label in record_points.keys():
        probe_id = f"v_{probe_label}"
        handle = sim.sample(
            0,
            probe_id,
            A.regular_schedule(dt_record_ms * A.units.ms),
        )
        probe_handles[probe_label] = handle
    logger.info("Probes recording at dt=%f ms", dt_record_ms)

    sim.record(A.spike_recording.all)
    logger.info(
        "Running simulation: T=%f ms, dt=%f ms...",
        self._parameters["T_ms"],
        self._parameters["dt_sim_ms"],
    )

    sim.run(
        self._parameters["T_ms"] * A.units.ms,
        self._parameters["dt_sim_ms"] * A.units.ms,
    )
    logger.info("Simulation completed")

    logger.info("Collecting voltage data from probes...")
    voltage_data_raw = {}
    for probe_label, handle in probe_handles.items():
        voltage_data_raw[probe_label] = sim.samples(handle)
    logger.info("Collected data from %d probes", len(voltage_data_raw))

    data_dict = {}
    probe_times = None
    for probe_label, samples in voltage_data_raw.items():
        if len(samples) > 0:
            data_array, location = samples[0]
            if probe_times is None:
                probe_times = data_array[:, 0]
            data_dict[probe_label] = data_array[:, 1]

    if probe_times is not None and len(data_dict) > 0:
        voltage_tsdframe = nap.TsdFrame(
            t=probe_times,
            d=np.column_stack([data_dict[k] for k in sorted(data_dict.keys())]),
            time_units="ms",
            columns=sorted(data_dict.keys()),
        )
    else:
        voltage_tsdframe = nap.TsdFrame(
            t=[],
            d=np.array([]).reshape(0, len(record_points)),
            time_units="ms",
            columns=list(record_points.keys()),
        )

    logger.info("Created voltage TsdFrame with shape %s", voltage_tsdframe.shape)

    simulation_results_dict = {
        "voltage_traces": voltage_tsdframe,
        "input_events": events,
        "synapses": self._synapses,
        "gap_junctions": self._gap_junctions,
        "record_points": record_points,
        "cell": self._build_cell_results["cell"],
        "morphology": self._build_cell_results["morphology"],
        "segment_tree": self._build_cell_results["segment_tree"],
        "decor": self._build_cell_results["decor"],
        "labels": self._build_cell_results["labels"],
        "cvp": self._build_cell_results["cvp"],
        "swc_filepath": self._swc_filepath,
        "synpts_filepath": self._synpts_filepath,
    }

    logger.info("Simulation complete")
    return SimulationResults(simulation_results_dict)

make_default_parameter_bank()

Returns a ParameterBank with default values. Call bank.sample() to obtain a ParameterSet for simulation.

Arbor ion channel defaults: Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM

Source code in toric_spines_sim/simulation/parameters.py

def make_default_parameter_bank() -> ParameterBank:
    """
    Returns a ParameterBank with default values.
    Call ``bank.sample()`` to obtain a ParameterSet for simulation.

    Arbor ion channel defaults:
    Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM
    K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM
    Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM


    """

    parameter_bank = ParameterBank(
        {
            "seed": IndependentScalarParameter(0),
            "T_ms": IndependentScalarParameter(1000.0),
            "delay_ms": IndependentScalarParameter(100.0),
            "discretization_um": IndependentScalarParameter(0.1),
            "dt_sim_ms": IndependentScalarParameter(0.02),
            "dt_record_ms": IndependentScalarParameter(0.1),
            "spine_tag": IndependentScalarParameter(3),
            "sink_tag": IndependentScalarParameter(5),
            "sink_tip_tag": IndependentScalarParameter(6),
            "sink_radii_scale": IndependentScalarParameter(1.0),
            "neck_radius_scale": IndependentScalarParameter(1.0),
            "temp_K": IndependentScalarParameter(
                280.0,
                is_sampled=False,
                range=(270.0, 320.0),
            ),
            "Vrest_mV": IndependentScalarParameter(
                -65.0,
                is_sampled=False,
                range=(-70.0, -60.0),
            ),
            "cm_uF_per_cm2": IndependentScalarParameter(
                1.0,
                is_sampled=False,
                range=(1.0, 20.0),
            ),
            "rL_ohm_cm": IndependentScalarParameter(
                35.0,
                is_sampled=False,
                range=(30.0, 200.0),
            ),
            "gj_weight": IndependentScalarParameter(1.0),
            "pas_leak_g_S_per_cm2": IndependentScalarParameter(
                0.001, is_sampled=False, range=(0.00001, 0.1),
            ),
            "pas_leak_e_mV": IndependentScalarParameter(-65.0),
            "tau_m_ms": DerivedScalarParameter(get_tau_m),
            "hh_leak_g_S_per_cm2": IndependentScalarParameter(
                0.0003, is_sampled=False, range=(0.00001, 0.01),
                ),
            "hh_leak_e_mV": IndependentScalarParameter(-54.3),
            "hh_tags": IndependentVectorParameter([]),
            "hh_scale": IndependentScalarParameter(1.0),
            "iclamp_locations": IndependentVectorParameter([]),
            "iclamp_amplitudes": IndependentVectorParameter([]),
            "iclamp_durations": IndependentVectorParameter([]),
            "K_revpot_mV": IndependentScalarParameter(-77.0),
            "K_intcon_mM": IndependentScalarParameter(54.4),
            "K_extcon_mM": IndependentScalarParameter(2.5),
            "K_gbar_S_per_cm2": IndependentScalarParameter(0.036),
            "Na_revpot_mV": IndependentScalarParameter(50.0),
            "Na_intcon_mM": IndependentScalarParameter(10.0),
            "Na_extcon_mM": IndependentScalarParameter(140.0),
            "Na_gbar_S_per_cm2": IndependentScalarParameter(0.12),
            "Ca_revpot_mV": IndependentScalarParameter(132.458),
            "Ca_intcon_mM": IndependentScalarParameter(0.00005),
            "Ca_extcon_mM": IndependentScalarParameter(2.0),
            "Ca_gbar_S_per_cm2": IndependentScalarParameter(0.0002),
            "ampa_gmax_uS": IndependentScalarParameter(
                0.002,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "ampa_tau_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 4.0),
            ),
            "ampa_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "nmda_gmax_uS": IndependentScalarParameter(
                0.002,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "nmda_tau_r_ms": IndependentScalarParameter(
                5.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "nmda_tau_d_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "nmda_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "gabaa_gmax_uS": IndependentScalarParameter(
                0.00008,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "gabaa_tau_ms": IndependentScalarParameter(
                10.0,
                is_sampled=False,
                range=(0.1, 50.0),
            ),
            "gabaa_e_mV": IndependentScalarParameter(
                -75.0,
                is_sampled=False,
                range=(-95.0, -60.0),
            ),
            "gabab_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "gabab_e_mV": IndependentScalarParameter(
                -95.0,
                is_sampled=False,
                range=(-100.0, -60.0),
            ),
            "gabab_tau_r_ms": IndependentScalarParameter(
                30.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "gabab_tau_d_ms": IndependentScalarParameter(
                200.0,
                is_sampled=False,
                range=(0.1, 500.0),
            ),
            "effexc_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "effexc_nmda_ratio": IndependentScalarParameter(
                0.50,
                is_sampled=False,
                range=(0.0, 1.0),
            ),
            "effexc_tau_ampa_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "effexc_tau_nmda_rise_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 20.0),
            ),
            "effexc_tau_nmda_decay_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 200.0),
            ),
            "effexc_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "effexc_mg_mM": IndependentScalarParameter(
                1.0,
                is_sampled=False,
                range=(0.0, 5.0),
            ),
        },
    )
    return parameter_bank

make_icx_parameter_bank_invitro()

ICx bank for comparison to Sanculi et al. patch-clamp (in vitro).

• Longitudinal resistivity (ohm cm): ~100 ohm cm, typical values 30 - 200
• Membrane capacitance (uF/cm^2): Cm = 8.559 +/- 50.475 (Sanculi; no in-vivo Cm)
• AMPA: gmax = 0.2-2.0 nS, tau = 2 ms, revpot = 0 mV
• NMDA: gmax = 0.2-2.0 nS, tau_r = 5 ms, tau_d = 50 ms, revpot = 0 mV
• GABA_A: gmax = 0.5-2.0 nS, tau = 2 ms, revpot = -70 mV

Arbor ion channel defaults: Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM

Source code in toric_spines_sim/simulation/parameters.py

def make_icx_parameter_bank_invitro() -> ParameterBank:
    """ICx bank for comparison to Sanculi et al. patch-clamp (in vitro).

    - Longitudinal resistivity (ohm cm): ~100 ohm cm, typical values 30 - 200
    - Membrane capacitance (uF/cm^2): Cm = 8.559 +/- 50.475 (Sanculi; no in-vivo Cm)
    - AMPA: gmax = 0.2-2.0 nS, tau = 2 ms, revpot = 0 mV
    - NMDA: gmax = 0.2-2.0 nS, tau_r = 5 ms, tau_d = 50 ms, revpot = 0 mV
    - GABA_A: gmax = 0.5-2.0 nS, tau = 2 ms, revpot = -70 mV

    Arbor ion channel defaults:
    Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM
    K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM
    Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM
    """

    parameter_bank = ParameterBank(
        {
            "seed": IndependentScalarParameter(0),
            "T_ms": IndependentScalarParameter(1000.0),
            "delay_ms": IndependentScalarParameter(100.0),
            "discretization_um": IndependentScalarParameter(0.1),
            "dt_sim_ms": IndependentScalarParameter(0.02),
            "dt_record_ms": IndependentScalarParameter(0.1),
            "spine_tag": IndependentScalarParameter(3),
            "sink_tag": IndependentScalarParameter(5),
            "sink_tip_tag": IndependentScalarParameter(6),
            "sink_radii_scale": IndependentScalarParameter(1.0),
            "neck_radius_scale": IndependentScalarParameter(1.0),
            "temp_K": IndependentScalarParameter(
                297.0,
                is_sampled=False,
                range=(270.0, 320.0),
            ),
            "Vrest_mV": IndependentScalarParameter(
                -67.6,
                is_sampled=False,
                range=(-80.0, -50.0),
            ),
            "cm_uF_per_cm2": IndependentScalarParameter(
                8.559,
                is_sampled=False,
                range=(1.0, 20.0),
            ),
            "rL_ohm_cm": IndependentScalarParameter(
                100.0,
                is_sampled=False,
                range=(30.0, 200.0),
            ),
            "gj_weight": IndependentScalarParameter(1.0),
            "pas_leak_g_S_per_cm2": IndependentScalarParameter(
                0.000144, is_sampled=False, range=(0.00001, 0.1),
            ),
            "pas_leak_e_mV": IndependentScalarParameter(-67.6),
            "tau_m_ms": DerivedScalarParameter(get_tau_m),
            "hh_leak_g_S_per_cm2": IndependentScalarParameter(
                0.0003, is_sampled=False, range=(0.00001, 0.01),
                ),
            "hh_leak_e_mV": IndependentScalarParameter(-54.3),
            "hh_tags": IndependentVectorParameter([]),
            "hh_scale": IndependentScalarParameter(1.0),
            "iclamp_locations": IndependentVectorParameter([]),
            "iclamp_amplitudes": IndependentVectorParameter([]),
            "iclamp_durations": IndependentVectorParameter([]),
            "K_revpot_mV": IndependentScalarParameter(-77.0),
            "K_intcon_mM": IndependentScalarParameter(54.4),
            "K_extcon_mM": IndependentScalarParameter(2.5),
            "K_gbar_S_per_cm2": IndependentScalarParameter(0.036),
            "Na_revpot_mV": IndependentScalarParameter(50.0),
            "Na_intcon_mM": IndependentScalarParameter(10.0),
            "Na_extcon_mM": IndependentScalarParameter(140.0),
            "Na_gbar_S_per_cm2": IndependentScalarParameter(0.12),
            "Ca_revpot_mV": IndependentScalarParameter(132.458),
            "Ca_intcon_mM": IndependentScalarParameter(0.00005),
            "Ca_extcon_mM": IndependentScalarParameter(2.0),
            "Ca_gbar_S_per_cm2": IndependentScalarParameter(0.0002),
            "ampa_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0002, 0.002),
            ),
            "ampa_tau_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 4.0),
            ),
            "ampa_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "nmda_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0002, 0.002),
            ),
            "nmda_tau_r_ms": IndependentScalarParameter(
                5.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "nmda_tau_d_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "nmda_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "gabaa_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0005, 0.002),
            ),
            "gabaa_tau_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 50.0),
            ),
            "gabaa_e_mV": IndependentScalarParameter(
                -70.0,
                is_sampled=False,
                range=(-95.0, -60.0),
            ),
            "gabab_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "gabab_e_mV": IndependentScalarParameter(
                -95.0,
                is_sampled=False,
                range=(-100.0, -60.0),
            ),
            "gabab_tau_r_ms": IndependentScalarParameter(
                30.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "gabab_tau_d_ms": IndependentScalarParameter(
                200.0,
                is_sampled=False,
                range=(0.1, 500.0),
            ),
            "effexc_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "effexc_nmda_ratio": IndependentScalarParameter(
                0.50,
                is_sampled=False,
                range=(0.0, 1.0),
            ),
            "effexc_tau_ampa_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "effexc_tau_nmda_rise_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 20.0),
            ),
            "effexc_tau_nmda_decay_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 200.0),
            ),
            "effexc_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "effexc_mg_mM": IndependentScalarParameter(
                1.0,
                is_sampled=False,
                range=(0.0, 5.0),
            ),
        },
    )
    return parameter_bank

make_icx_parameter_bank_invivo()

ICx bank for comparison to living-owl intracellular recordings. Peña and Konishi (2002, Journal of Neuroscience): mean rest potential = -67.6 ± 9.3 mV (n = 75 ICx neurons, sharp electrode). Body temperature = ~313 K. Input resistance is lower in vivo than in slice (~3×); leak = Sanculi g_pas × 3 = 0.00042 S/cm². Cm is still the Sanculi value (no in-vivo capacitance).

• Longitudinal resistivity (ohm cm): ~100 ohm cm, typical values 30 - 200
• Membrane capacitance (uF/cm^2): Cm = 8.559 +/- 50.475 (Sanculi; no in-vivo Cm)
• AMPA: gmax = 0.2-2.0 nS, tau = 2 ms, revpot = 0 mV
• NMDA: gmax = 0.2-2.0 nS, tau_r = 5 ms, tau_d = 50 ms, revpot = 0 mV
• GABA_A: gmax = 0.5-2.0 nS, tau = 2 ms, revpot = -70 mV

Arbor ion channel defaults: Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM

Source code in toric_spines_sim/simulation/parameters.py

def make_icx_parameter_bank_invivo() -> ParameterBank:
    """
    ICx bank for comparison to living-owl intracellular recordings.
    Peña and Konishi (2002, Journal of Neuroscience): 
    mean rest potential = -67.6 ± 9.3 mV (n = 75 ICx neurons, sharp electrode). 
    Body temperature = ~313 K. 
    Input resistance is lower in vivo than in slice (~3×); 
    leak = Sanculi g_pas × 3 = 0.00042 S/cm². 
    Cm is still the Sanculi value (no in-vivo capacitance).

    - Longitudinal resistivity (ohm cm): ~100 ohm cm, typical values 30 - 200
    - Membrane capacitance (uF/cm^2): Cm = 8.559 +/- 50.475 (Sanculi; no in-vivo Cm)
    - AMPA: gmax = 0.2-2.0 nS, tau = 2 ms, revpot = 0 mV
    - NMDA: gmax = 0.2-2.0 nS, tau_r = 5 ms, tau_d = 50 ms, revpot = 0 mV
    - GABA_A: gmax = 0.5-2.0 nS, tau = 2 ms, revpot = -70 mV

    Arbor ion channel defaults:
    Na: revpot = 50.0 mV, int_conc = 10.0 mM, ext_conc = 140.0 mM
    K: revpot = -77.0 mV, int_conc = 54.4 mM, ext_conc = 2.5 mM
    Ca: revpot = 132.458 mV, int_conc = 0.00005 mM, ext_conc = 2 mM

    """

    parameter_bank = ParameterBank(
        {
            "seed": IndependentScalarParameter(0),
            "T_ms": IndependentScalarParameter(1000.0),
            "delay_ms": IndependentScalarParameter(100.0),
            "discretization_um": IndependentScalarParameter(0.1),
            "dt_sim_ms": IndependentScalarParameter(0.02),
            "dt_record_ms": IndependentScalarParameter(0.1),
            "spine_tag": IndependentScalarParameter(3),
            "sink_tag": IndependentScalarParameter(5),
            "sink_tip_tag": IndependentScalarParameter(6),
            "sink_radii_scale": IndependentScalarParameter(1.0),
            "neck_radius_scale": IndependentScalarParameter(1.0),
            "temp_K": IndependentScalarParameter(
                313.0,
                is_sampled=False,
                range=(270.0, 320.0),
            ),
            "Vrest_mV": IndependentScalarParameter(
                -67.6,
                is_sampled=False,
                range=(-80.0, -50.0),
            ),
            "cm_uF_per_cm2": IndependentScalarParameter(
                8.559,
                is_sampled=False,
                range=(1.0, 20.0),
            ),
            "rL_ohm_cm": IndependentScalarParameter(
                100.0,
                is_sampled=False,
                range=(30.0, 200.0),
            ),
            "gj_weight": IndependentScalarParameter(1.0),
            "pas_leak_g_S_per_cm2": IndependentScalarParameter(
                0.00042, is_sampled=False, range=(0.00001, 0.1),
            ),
            "pas_leak_e_mV": IndependentScalarParameter(-67.6),
            "tau_m_ms": DerivedScalarParameter(get_tau_m),
            "hh_leak_g_S_per_cm2": IndependentScalarParameter(
                0.0003, is_sampled=False, range=(0.00001, 0.01),
                ),
            "hh_leak_e_mV": IndependentScalarParameter(-54.3),
            "hh_tags": IndependentVectorParameter([]),
            "hh_scale": IndependentScalarParameter(1.0),
            "iclamp_locations": IndependentVectorParameter([]),
            "iclamp_amplitudes": IndependentVectorParameter([]),
            "iclamp_durations": IndependentVectorParameter([]),
            "K_revpot_mV": IndependentScalarParameter(-77.0),
            "K_intcon_mM": IndependentScalarParameter(54.4),
            "K_extcon_mM": IndependentScalarParameter(2.5),
            "K_gbar_S_per_cm2": IndependentScalarParameter(0.036),
            "Na_revpot_mV": IndependentScalarParameter(50.0),
            "Na_intcon_mM": IndependentScalarParameter(10.0),
            "Na_extcon_mM": IndependentScalarParameter(140.0),
            "Na_gbar_S_per_cm2": IndependentScalarParameter(0.12),
            "Ca_revpot_mV": IndependentScalarParameter(132.458),
            "Ca_intcon_mM": IndependentScalarParameter(0.00005),
            "Ca_extcon_mM": IndependentScalarParameter(2.0),
            "Ca_gbar_S_per_cm2": IndependentScalarParameter(0.0002),
            "ampa_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0002, 0.002),
            ),
            "ampa_tau_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 4.0),
            ),
            "ampa_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "nmda_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0002, 0.002),
            ),
            "nmda_tau_r_ms": IndependentScalarParameter(
                5.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "nmda_tau_d_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "nmda_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "gabaa_gmax_uS": IndependentScalarParameter(
                0.001,
                is_sampled=False,
                range=(0.0005, 0.002),
            ),
            "gabaa_tau_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 50.0),
            ),
            "gabaa_e_mV": IndependentScalarParameter(
                -70.0,
                is_sampled=False,
                range=(-95.0, -60.0),
            ),
            "gabab_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "gabab_e_mV": IndependentScalarParameter(
                -95.0,
                is_sampled=False,
                range=(-100.0, -60.0),
            ),
            "gabab_tau_r_ms": IndependentScalarParameter(
                30.0,
                is_sampled=False,
                range=(0.1, 100.0),
            ),
            "gabab_tau_d_ms": IndependentScalarParameter(
                200.0,
                is_sampled=False,
                range=(0.1, 500.0),
            ),
            "effexc_gmax_uS": IndependentScalarParameter(
                0.00010,
                is_sampled=False,
                range=(0.0, 0.5),
            ),
            "effexc_nmda_ratio": IndependentScalarParameter(
                0.50,
                is_sampled=False,
                range=(0.0, 1.0),
            ),
            "effexc_tau_ampa_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 10.0),
            ),
            "effexc_tau_nmda_rise_ms": IndependentScalarParameter(
                2.0,
                is_sampled=False,
                range=(0.1, 20.0),
            ),
            "effexc_tau_nmda_decay_ms": IndependentScalarParameter(
                50.0,
                is_sampled=False,
                range=(0.1, 200.0),
            ),
            "effexc_e_mV": IndependentScalarParameter(
                0.0,
                is_sampled=False,
                range=(-70.0, 0.0),
            ),
            "effexc_mg_mM": IndependentScalarParameter(
                1.0,
                is_sampled=False,
                range=(0.0, 5.0),
            ),
        },
    )
    return parameter_bank

random_axon_events(n_synapses, n_axons, mod_freq_hz, peak_rate_hz=None, peak_rate_range_hz=None, phase_range_rad=(0.0, 2.0 * np.pi), seed=42)

Create random sine rate curves and a random synapse-to-axon partition.

Provide exactly one of peak_rate_hz or peak_rate_range_hz.

Parameters:

Name	Type	Description	Default
n_synapses	int	Total number of synapses.	required
n_axons	int	Number of axons (must be <= n_synapses).	required
mod_freq_hz	float	Common sine frequency (Hz).	required
peak_rate_hz	float	Peak rate for every axon (Hz).	None
peak_rate_range_hz	tuple of float	(min, max) from which each axon's peak is drawn.	None
phase_range_rad	tuple of float	Phase draw range in radians. Default (0, 2π).	(0.0, 2.0 * pi)
seed	int	RNG seed.	42

Returns:

Name	Type	Description
rate_curves	list of SineRateCurve
n_synapses_per_axon	list of int

Source code in toric_spines_sim/simulation/input.py

def random_axon_events(
    n_synapses: int,
    n_axons: int,
    mod_freq_hz: float,
    peak_rate_hz: float = None,
    peak_rate_range_hz: Tuple[float, float] = None,
    phase_range_rad: Tuple[float, float] = (0.0, 2.0 * np.pi),
    seed: int = 42,
) -> Tuple[List[SineRateCurve], List[int]]:
    """Create random sine rate curves and a random synapse-to-axon partition.

    Provide exactly one of ``peak_rate_hz`` or ``peak_rate_range_hz``.

    Parameters
    ----------
    n_synapses : int
        Total number of synapses.
    n_axons : int
        Number of axons (must be ``<= n_synapses``).
    mod_freq_hz : float
        Common sine frequency (Hz).
    peak_rate_hz : float, optional
        Peak rate for every axon (Hz).
    peak_rate_range_hz : tuple of float, optional
        ``(min, max)`` from which each axon's peak is drawn.
    phase_range_rad : tuple of float
        Phase draw range in radians. Default ``(0, 2π)``.
    seed : int
        RNG seed.

    Returns
    -------
    rate_curves : list of SineRateCurve
    n_synapses_per_axon : list of int
    """
    if n_axons > n_synapses:
        raise ValueError(f"n_axons ({n_axons}) must be <= n_synapses ({n_synapses})")

    rng = np.random.default_rng(seed)

    # Randomly assign synapses to axons
    # Start with at least 1 synapse per axon, then distribute remaining
    base_assignment = np.arange(n_axons) % n_synapses
    remaining = n_synapses - n_axons
    extra_assignments = rng.integers(0, n_axons, size=remaining)
    assignments = np.concatenate([base_assignment, extra_assignments])
    rng.shuffle(assignments)

    # Count synapses per axon
    n_synapses_per_axon = [int((assignments == i).sum()) for i in range(n_axons)]

    # Create random rate curves for each axon
    rate_curves = []
    for i in range(n_axons):
        if peak_rate_hz is not None:
            peak_rate = peak_rate_hz
        elif peak_rate_range_hz is not None:
            peak_rate = rng.uniform(peak_rate_range_hz[0], peak_rate_range_hz[1])
        else:
            raise ValueError("Either peak_rate_hz or peak_rate_range_hz must be provided")
        phase = rng.uniform(phase_range_rad[0], phase_range_rad[1])
        rate_curves.append(
            SineRateCurve(peak_rate_hz=peak_rate, freq_hz=mod_freq_hz, phase_rad=phase)
        )

    logger.info(
        f"Created {n_axons} axons for {n_synapses} synapses "
        f"(synapses_per_axon: {n_synapses_per_axon})"
    )
    return rate_curves, n_synapses_per_axon

load_axon_events_from_file(axon_assignment_file, axon_rates_hz)

Build per-axon flat rate curves from a synapse assignment file.

Each line is one axon: comma-separated 1-based synapse indices.

Parameters:

Name	Type	Description	Default
axon_assignment_file	path - like	Assignment file (one axon per line).	required
axon_rates_hz	list of float	Rate in Hz for each axon. Length must match the number of lines. Use 0.0 to silence an axon.	required

Returns:

Name	Type	Description
rate_curves	list of FlatRateCurve
n_synapses_per_axon	list of int
axon_synapses	list of list of int	0-based synapse indices per axon (for remap_axon_channel_events_to_synapses).

Source code in toric_spines_sim/simulation/input.py

def load_axon_events_from_file(
    axon_assignment_file: Path,
    axon_rates_hz: List[float],
) -> Tuple[List[FlatRateCurve], List[int], List[List[int]]]:
    """Build per-axon flat rate curves from a synapse assignment file.

    Each line is one axon: comma-separated **1-based** synapse indices.

    Parameters
    ----------
    axon_assignment_file : path-like
        Assignment file (one axon per line).
    axon_rates_hz : list of float
        Rate in Hz for each axon. Length must match the number of lines.
        Use ``0.0`` to silence an axon.

    Returns
    -------
    rate_curves : list of FlatRateCurve
    n_synapses_per_axon : list of int
    axon_synapses : list of list of int
        0-based synapse indices per axon (for ``remap_axon_channel_events_to_synapses``).
    """
    # Read synapse assignments from file
    axon_synapses = []
    with open(axon_assignment_file, 'r') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            # Parse comma-separated indices (1-based), convert to 0-based
            synapse_indices = [int(x.strip()) - 1 for x in line.split(',')]
            axon_synapses.append(synapse_indices)

    n_axons = len(axon_synapses)
    if len(axon_rates_hz) != n_axons:
        raise ValueError(
            f"axon_rates_hz length ({len(axon_rates_hz)}) must match "
            f"number of axons in file ({n_axons})"
        )

    # Create flat rate curves for each axon
    rate_curves = []
    for rate_hz in axon_rates_hz:
        rate_curves.append(FlatRateCurve(rate_hz=rate_hz))

    # Count synapses per axon
    n_synapses_per_axon = [len(synapses) for synapses in axon_synapses]

    total_synapses = sum(n_synapses_per_axon)
    logger.info(
        f"Loaded {n_axons} axons for {total_synapses} synapses from "
        f"{axon_assignment_file} (synapses_per_axon: {n_synapses_per_axon})"
    )

    return rate_curves, n_synapses_per_axon, axon_synapses

remap_axon_channel_events_to_synapses(events_tsgroup, axon_synapses, n_synapses=None)

Map axon-ordered event channels onto synapse point-file indices.

Event generators fan out channels in axon order (axon 0 synapses, then axon 1, …), but TSRecipe maps TsGroup index i to synapse syn_i. This scatters axon-channel events onto the correct indices.

Parameters:

Name	Type	Description	Default
events_tsgroup	TsGroup	Channels ordered by axon assignment (generator output).	required
axon_synapses	list of list of int	Per-axon 0-based synapse indices (from load_axon_events_from_file).	required
n_synapses	int	Total synapses. Defaults to max(index) + 1 over assignments.	None

Returns:

Type	Description
TsGroup	Indexed 0 .. n_synapses-1 with labels syn_0, syn_1, …

Examples:

>>> # Axon 0 hits synapses 2 then 0; axon 1 hits synapse 1
>>> axon_synapses = [[2, 0], [1]]
>>> remapped = remap_axon_channel_events_to_synapses(events, axon_synapses)
>>> remapped.get_info("label")[0]
'syn_0'

Source code in toric_spines_sim/simulation/input.py

def remap_axon_channel_events_to_synapses(
    events_tsgroup: nap.TsGroup,
    axon_synapses: List[List[int]],
    n_synapses: Optional[int] = None,
) -> nap.TsGroup:
    """Map axon-ordered event channels onto synapse point-file indices.

    Event generators fan out channels in axon order (axon 0 synapses, then
    axon 1, …), but ``TSRecipe`` maps TsGroup index ``i`` to synapse
    ``syn_i``. This scatters axon-channel events onto the correct indices.

    Parameters
    ----------
    events_tsgroup : pynapple.TsGroup
        Channels ordered by axon assignment (generator output).
    axon_synapses : list of list of int
        Per-axon 0-based synapse indices (from ``load_axon_events_from_file``).
    n_synapses : int, optional
        Total synapses. Defaults to ``max(index) + 1`` over assignments.

    Returns
    -------
    pynapple.TsGroup
        Indexed ``0 .. n_synapses-1`` with labels ``syn_0``, ``syn_1``, …

    Examples
    --------
    >>> # Axon 0 hits synapses 2 then 0; axon 1 hits synapse 1
    >>> axon_synapses = [[2, 0], [1]]
    >>> remapped = remap_axon_channel_events_to_synapses(events, axon_synapses)  # doctest: +SKIP
    >>> remapped.get_info("label")[0]
    'syn_0'
    """
    expected_channels = sum(len(synapses) for synapses in axon_synapses)
    if len(events_tsgroup) != expected_channels:
        raise ValueError(
            f"events_tsgroup has {len(events_tsgroup)} channels, expected "
            f"{expected_channels} from axon_synapses"
        )

    if n_synapses is None:
        if not any(axon_synapses):
            raise ValueError("axon_synapses is empty; cannot infer n_synapses")
        n_synapses = max(syn_idx for synapses in axon_synapses for syn_idx in synapses) + 1

    remapped_by_index: dict[int, nap.Ts] = {}
    channel_idx = 0
    for synapse_indices in axon_synapses:
        for syn_idx in synapse_indices:
            if syn_idx in remapped_by_index:
                raise ValueError(
                    f"Synapse index {syn_idx} appears in multiple axon assignments"
                )
            if syn_idx < 0 or syn_idx >= n_synapses:
                raise ValueError(
                    f"Synapse index {syn_idx} out of range for n_synapses={n_synapses}"
                )
            remapped_by_index[syn_idx] = events_tsgroup[channel_idx]
            channel_idx += 1

    # Build 0..n-1 in order so TsGroup label metadata aligns with synapse indices.
    remapped = {
        syn_idx: remapped_by_index.get(
            syn_idx, nap.Ts(t=[], time_units="ms")
        )
        for syn_idx in range(n_synapses)
    }

    labels = [f"syn_{i}" for i in range(n_synapses)]
    return nap.TsGroup(
        remapped,
        label=labels,
        time_support=events_tsgroup.time_support,
    )