model

Cell morphology and recipe definitions for toric spine models.

toric_spines_sim.model

Cell morphology and recipe definitions for toric spine models.

TSModel(swc_path, synapses, gap_junctions, record_points, parameters) dataclass

Morphology specification for toric spine models.

Attributes:

Name	Type	Description
swc_path	Path	Path to the base SWC morphology for the toric spine.
synapses	dict[str, SynapsePoint]	Dictionary of synapse labels to instances of SynapsePoint.
gap_junctions	dict[str, GapJunctionPoint]	Dictionary of gap junction labels to instances of GapJunctionPoint.
record_points	dict[str, Tuple[float, float, float]]	Dictionary of record site labels to 3D points (x,y,z).
parameters	ParameterSet	Sampled simulation parameters (from ParameterBank.sample()).

build_cell()

Construct an Arbor cable cell from the SWC and explicit placements.

Returns:

Type	Description
dict	Intermediate objects keyed by: cell : arbor.cable_cell morphology : arbor.morphology segment_tree : arbor.segment_tree decor : arbor.decor labels : arbor.label_dict cvp : CV policy custom_catalogue : loaded NMODL catalogue

Examples:

>>> build_result = model.build_cell()
>>> cell = build_result["cell"]

Source code in toric_spines_sim/model/model.py

def build_cell(self) -> dict[str, object]:
    """Construct an Arbor cable cell from the SWC and explicit placements.

    Returns
    -------
    dict
        Intermediate objects keyed by:

        - ``cell`` : ``arbor.cable_cell``
        - ``morphology`` : ``arbor.morphology``
        - ``segment_tree`` : ``arbor.segment_tree``
        - ``decor`` : ``arbor.decor``
        - ``labels`` : ``arbor.label_dict``
        - ``cvp`` : CV policy
        - ``custom_catalogue`` : loaded NMODL catalogue

    Examples
    --------
    >>> build_result = model.build_cell()  # doctest: +SKIP
    >>> cell = build_result["cell"]
    """
    # Load pre-built custom mechanism catalogue (ampasyn, nmdasyn, hhnotemp, ...)
    custom_catalogue = check_catalogue()

    # Load swc morphology (fallback for older Arbor without raw=True)
    logger.info(
        "Building cable_cell from %s with discretization %.3f um",
        self.swc_path,
        self.parameters["discretization_um"],
    )
    loaded_morphology = A.load_swc_arbor(str(self.swc_path))
    original_segment_tree = getattr(
        loaded_morphology, "segment_tree", loaded_morphology
    )
    segment_tree = original_segment_tree
    if self.parameters["sink_radii_scale"] != 1.0:
        segment_tree = scale_radii_in_segment_tree_by_tag(
            segment_tree,
            self.parameters["sink_radii_scale"],
            self.parameters["sink_tag"],
        )
        logger.debug(
            "Scaled sink radii by %.4f",
            self.parameters["sink_radii_scale"],
        )
    if self.parameters["neck_radius_scale"] != 1.0:
        neck_xyz = neck_point_from_swc_file(self.swc_path)
        segment_tree = scale_one_radius_in_segment_tree_by_coordinates(
            segment_tree,
            self.parameters["neck_radius_scale"],
            neck_xyz,
        )
        logger.debug(
            "Scaled neck radius at %s by %.4f",
            neck_xyz,
            self.parameters["neck_radius_scale"],
        )

    morphology = A.morphology(segment_tree)

    hh_tags = [int(tag) for tag in self.parameters["hh_tags"]]
    hh_scale = (
        float(self.parameters["hh_scale"])
        if "hh_scale" in self.parameters.index
        else 0.0
    )
    label_map = {
        "all": "(all)",
        "root": "(root)",
        "sink": f"(tag {self.parameters['sink_tag']})",
        "spine": f"(tag {self.parameters['spine_tag']})",
    }
    if len(hh_tags) > 0 and hh_scale > 0:
        label_map["hh_area"] = join_tags_dsl(hh_tags)

    decor = A.decor()

    # DENSITY MECHANISMS ===============================
    # Passive leak everywhere (independent of HH leak parameters).
    pas_mechanism_name = f"pas/e={self.parameters['pas_leak_e_mV']}"
    decor.paint(
        label_map["all"],
        A.density(
            pas_mechanism_name,
            {
                "g": self.parameters["pas_leak_g_S_per_cm2"],
            },
        ),
    )
    # Hodgkin-Huxley (only if hh_tags are specified and hh_scale > 0)
    use_hh_notemp = True
    if "hh_area" in label_map:
        hh_mechanism_name = "hhnotemp" if use_hh_notemp else "hh"
        decor.paint(
            label_map["hh_area"],
            A.density(
                hh_mechanism_name,
                {
                    "gnabar": self.parameters["Na_gbar_S_per_cm2"] * hh_scale,
                    "gkbar": self.parameters["K_gbar_S_per_cm2"] * hh_scale,
                    "gl": self.parameters["hh_leak_g_S_per_cm2"] * hh_scale,
                    "el": self.parameters["hh_leak_e_mV"],
                },
            ),
        )
        logger.info(
            "Painted %s on tags %s with scale %.4g "
            "(gnabar=%.6g, gkbar=%.6g g/cm2)",
            hh_mechanism_name,
            hh_tags,
            hh_scale,
            self.parameters["Na_gbar_S_per_cm2"] * hh_scale,
            self.parameters["K_gbar_S_per_cm2"] * hh_scale,
        )

    # Map synapse 3D points to closest Arbor locations; place individually under unique labels.
    piecewise_placer = A.place_pwlin(morphology)
    if self.synapses:
        logger.debug("Placing %d synapses", len(self.synapses))
        for syn_label, synapse in self.synapses.items():
            x, y, z = synapse.location
            location, _ = piecewise_placer.closest(float(x), float(y), float(z))
            location_expr = f"(location {location.branch} {location.pos:.6f})"

            decor.place(
                location_expr,
                A.synapse(synapse.mechanism, synapse.mechanism_params),
                syn_label,
            )
            label_map[syn_label] = location_expr
            logger.debug(
                "Placed synapse %s at %s (mechanism=%s, params=%s)",
                syn_label,
                location_expr,
                synapse.mechanism,
                synapse.mechanism_params,
            )
    else:
        logger.info("No synapses provided.")

    # Restore cycles / extra necks: place two junction labels on the
    # distinct SWC samples named by index_pair (not a shared closest XYZ).
    if self.gap_junctions:
        node_ids = []
        for gj in self.gap_junctions.values():
            node_ids.extend(gj.index_pair)
        node_locations = arbor_locations_for_swc_nodes(
            self.swc_path, morphology, segment_tree, node_ids
        )
        logger.debug("Placing %d gap junctions", len(self.gap_junctions))
        for gj_label, gj in self.gap_junctions.items():
            node_i, node_j = gj.index_pair
            gj_label_a = f"{gj_label}_a"
            gj_label_b = f"{gj_label}_b"
            location_i = node_locations[node_i]
            location_j = node_locations[node_j]
            location_expr_i = f"(location {location_i.branch} {location_i.pos:.6f})"
            location_expr_j = f"(location {location_j.branch} {location_j.pos:.6f})"
            if (
                location_i.branch == location_j.branch
                and abs(location_i.pos - location_j.pos) < 1e-9
            ):
                logger.warning(
                    "Gap junction %s maps both SWC nodes %s and %s to %s; "
                    "the connection may be electrically inert",
                    gj_label,
                    node_i,
                    node_j,
                    location_expr_i,
                )
            decor.place(location_expr_i, A.junction("gj"), gj_label_a)
            decor.place(location_expr_j, A.junction("gj"), gj_label_b)
            label_map[gj_label_a] = location_expr_i
            label_map[gj_label_b] = location_expr_j
            label_map[gj_label] = location_expr_i
            logger.debug(
                "Placed gap junction %s at %s / %s (SWC nodes %s, %s)",
                gj_label,
                location_expr_i,
                location_expr_j,
                node_i,
                node_j,
            )
    else:
        logger.info("No gap junctions provided.")

    if self.record_points:
        logger.debug("Placing %d record points", len(self.record_points))
        for record_label, record_point in self.record_points.items():
            x, y, z = record_point
            location, _ = piecewise_placer.closest(float(x), float(y), float(z))
            location_expr = f"(location {location.branch} {location.pos:.6f})"
            label_map[record_label] = location_expr
            logger.debug(
                "Placed record point %s at %s", record_label, location_expr
            )
    else:
        logger.info("No record points provided.")

    labels = A.label_dict(label_map)
    cvp = A.cv_policy_max_extent(
        self.parameters["discretization_um"] * U.um
    )
    cell = A.cable_cell(morphology, decor, labels, cvp)
    logger.info("Built cable_cell with %d labels", len(label_map))
    output = {
        "cell": cell,
        "morphology": morphology,
        "segment_tree": segment_tree,
        "decor": decor,
        "labels": labels,
        "cvp": cvp,
        "custom_catalogue": custom_catalogue,
    }
    return output

TSRecipe(cell, synapses, gap_junctions, record_points, events=None, parameters=None, custom_catalogue=None)

Bases: recipe

Arbor recipe for a single toric-spine cell with per-synapse events.

Parameters:

Name	Type	Description	Default
cell	cable_cell	Cell from TSModel.build_cell()["cell"].	required
synapses	dict[str, SynapsePoint]	Synapse labels used as Arbor place tags (typically syn_0, …).	required
gap_junctions	dict[str, GapJunctionPoint]	Gap junctions for CYCLE_BREAK / MULTI_NECK reconnects.	required
record_points	dict[str, tuple]	Voltage probe labels and XYZ positions.	required
events	TsGroup or dict[str, list[float]]	Event times in milliseconds. A TsGroup is mapped by index to list(synapses.keys()) order, not by channel label. A dict maps synapse labels to time lists directly.	None
parameters	ParameterSet	Required sampled bank (ions, leak, temperature, …).	None
custom_catalogue	catalogue	Extra NMODL mechanisms (ampa, nmda, hhnotemp, …).	None

Source code in toric_spines_sim/model/recipe.py

def __init__(
    self,
    cell: A.cable_cell,
    synapses: Dict[str, SynapsePoint],
    gap_junctions: Dict[str, GapJunctionPoint],
    record_points: Dict[str, Tuple[float, float, float]],
    events: Optional[Union[nap.TsGroup, Dict[str, List[float]]]] = None,
    parameters: Optional[ParameterSet] = None,
    custom_catalogue: Optional[A.catalogue] = None,
):
    super().__init__()
    if parameters is None:
        raise TypeError(
            "TSRecipe requires a ParameterSet (ion concentrations and "
            "reversal potentials). Pass parameters= from a sampled parameter bank."
        )
    self._cell = cell
    self._synapses = synapses
    self._gap_junctions = gap_junctions
    self._record_points = record_points
    self._parameters = parameters
    self._custom_catalogue = custom_catalogue

    # Convert TsGroup to dict if needed
    if isinstance(events, nap.TsGroup):
        self._events_ms = self._tsgroup_to_dict(events, list(synapses.keys()))
    else:
        self._events_ms = events

    Vrest = self._parameters["Vrest_mV"] * U.mV
    tempK = self._parameters["temp_K"] * U.Kelvin
    cm = self._parameters["cm_uF_per_cm2"] * U.uF / U.cm2
    rL = self._parameters["rL_ohm_cm"] * U.Ohm * U.cm

    # Global properties (passive defaults; catalog)
    self._gprop = A.cable_global_properties()
    self._gprop.catalogue = A.default_catalogue()

    # Extend with custom catalogue if provided 
    if self._custom_catalogue is not None:
        logger.info("Extending catalogue with custom mechanisms")
        self._gprop.catalogue.extend(self._custom_catalogue, "")

    self._gprop.set_property(Vm=Vrest, cm=cm, rL=rL, tempK=tempK)
    self._gprop.set_ion(
        "ca",
        valence=2,
        int_con=self._parameters["Ca_intcon_mM"] * U.mM,
        ext_con=self._parameters["Ca_extcon_mM"] * U.mM,
        rev_pot=self._parameters["Ca_revpot_mV"] * U.mV,
    )
    self._gprop.set_ion(
        "na",
        valence=1,
        int_con=self._parameters["Na_intcon_mM"] * U.mM,
        ext_con=self._parameters["Na_extcon_mM"] * U.mM,
        rev_pot=self._parameters["Na_revpot_mV"] * U.mV,
    )
    self._gprop.set_ion(
        "k",
        valence=1,
        int_con=self._parameters["K_intcon_mM"] * U.mM,
        ext_con=self._parameters["K_extcon_mM"] * U.mM,
        rev_pot=self._parameters["K_revpot_mV"] * U.mV,
    )
    logger.info(
        "Initialized TSRecipe with %d synapses, %d gap junctions, %d record points",
        len(self._synapses) if self._synapses else 0,
        len(self._gap_junctions) if self._gap_junctions else 0,
        len(self._record_points) if self._record_points else 0,
    )

    # Precompute explicit schedules per-synapse if provided
    self._evgens: List[A.event_generator] = []
    if self._events_ms is not None and self._synapses:
        syn_labels = list(self._synapses.keys())
        if len(self._events_ms) != len(syn_labels):
            logger.error(
                "events length %d must match number of synapses %d",
                len(self._events_ms),
                len(syn_labels),
            )
            raise ValueError(
                f"events length {len(self._events_ms)} must match number of synapses {len(syn_labels)}"
            )
        for syn_label, time_list in self._events_ms.items():
            times_quantities = [float(t) * U.ms for t in (time_list or [])]
            schedule = A.explicit_schedule(times_quantities)
            self._evgens.append(A.event_generator(syn_label, 1.0, schedule))  # gmax is set in the synapse mechanism parameters
        logger.debug("Created %d event generators", len(self._evgens))

GapJunctionPoint(index_pair, location, weight, location_b=None) dataclass

One electrical reconnect between two SWC samples.

Attributes:

Name	Type	Description
index_pair	Tuple[int, int]	SWC node IDs (i, j) to join.
location	Tuple[float, float, float]	XYZ of node i (µm or px, matching the SWC). Display only; electrical endpoints are index_pair.
weight	float	Arbor gap-junction weight (dimensionless conductance scale).
location_b	Optional[Tuple[float, float, float]]	XYZ of node j when it differs from i; otherwise None. Display only: TSModel places the junction from index_pair, not from these coordinates.

SynapsePoint(location, model, mechanism, synapse_params, mechanism_params) dataclass

One placed synapse: XYZ, mechanism name, and parameter dicts.

synapse_params uses human-readable keys (gmax_uS, tau_ms, …). mechanism_params is the same values renamed for Arbor (gmax, tau, …). TSModel applies mechanism_params.

Attributes:

Name	Type	Description
location	tuple of float	XYZ in the same units as the SWC.
model	str	Registry key (ampa, nmda, gabaa, gabab, effexc).
mechanism	str	NMODL mechanism name (e.g. ampasyn).
synapse_params, mechanism_params	dict	Conductance/time-constant values; see above.

Examples:

>>> synapse_population = SynapsePopulation.from_file("TS1_synpts.txt", "ampa", params)
>>> synapse_population.synapses["syn_0"].mechanism_params["gmax"]

SynapsePopulation(model, locations, global_parameters, parameter_override=None, label_prefix='syn')

A homogeneous population of synapses sharing one model type.

Each synapse is assigned parameters by sampling an internal ParameterBank built from global_parameters. By default all parameters have is_sampled=False, so each synapse receives the same values.

Examples:

>>> synapse_population = SynapsePopulation.from_file("TS1_synpts.txt", "ampa", params)
>>> synapse_population.synapses["syn_0"].mechanism
'ampasyn'

Source code in toric_spines_sim/model/synapse.py

def __init__(
    self,
    model: str,
    locations: Sequence[Location],
    global_parameters: ParameterSet,
    parameter_override: Optional[ParameterBank] = None,
    label_prefix: str = "syn",
):
    spec = _validate_model(model)
    self.model = model
    self.mechanism: str = spec["mechanism"]  # type: ignore[assignment]
    self._parameter_bank = _build_parameter_bank(
        model, global_parameters, parameter_override
    )
    self.synapses: Dict[str, SynapsePoint] = {}

    for index, location in enumerate(locations):
        sampled = self._parameter_bank.sample()
        synapse_params = _extract_synapse_params(sampled, model)
        mechanism_params = _build_mechanism_params(synapse_params, model)
        label = f"{label_prefix}_{index}"
        self.synapses[label] = SynapsePoint(
            location=location,
            model=model,
            mechanism=self.mechanism,
            synapse_params=synapse_params,
            mechanism_params=mechanism_params,
        )

    logger.debug(
        "Built %d %s synapses (mechanism=%s)",
        len(self.synapses),
        model,
        self.mechanism,
    )

from_file(points_file, model, global_parameters, parameter_override=None, label_prefix='syn') classmethod

Load XYZ locations from a file and build a synapse population.

Source code in toric_spines_sim/model/synapse.py

@classmethod
def from_file(
    cls,
    points_file: Union[str, Path],
    model: str,
    global_parameters: ParameterSet,
    parameter_override: Optional[ParameterBank] = None,
    label_prefix: str = "syn",
) -> SynapsePopulation:
    """Load XYZ locations from a file and build a synapse population."""
    locations = _load_locations(points_file)
    population = cls(
        model=model,
        locations=locations,
        global_parameters=global_parameters,
        parameter_override=parameter_override,
        label_prefix=label_prefix,
    )
    logger.debug(
        "Prepared %d synapses from %s (model=%s)",
        len(population.synapses),
        points_file,
        model,
    )
    return population

merge(*populations) staticmethod

Combine multiple populations into one label -> SynapsePoint dict.

Source code in toric_spines_sim/model/synapse.py

@staticmethod
def merge(*populations: SynapsePopulation) -> Dict[str, SynapsePoint]:
    """Combine multiple populations into one label -> SynapsePoint dict."""
    merged: Dict[str, SynapsePoint] = {}
    for population in populations:
        for label, synapse in population.synapses.items():
            if label in merged:
                raise ValueError(f"Duplicate synapse label: {label}")
            merged[label] = synapse
    return merged

check_catalogue(path=None)

Load the custom NMODL catalogue, or raise if it is missing or unusable.

Checks that custom-catalogue.so exists, that Arbor can load it, and that every mechanism in mechanisms/my_catalogue/*.mod is present. Rebuild after changing .mod files, upgrading Arbor, or changing OS/compiler; do not copy a .so between machines.

Parameters:

Name	Type	Description	Default
path	Path	Catalogue file. Default: CUSTOM_CATALOGUE_PATH.	None

Returns:

Type	Description
catalogue

Source code in toric_spines_sim/model/model.py

def check_catalogue(path: Path | None = None) -> A.catalogue:
    """Load the custom NMODL catalogue, or raise if it is missing or unusable.

    Checks that ``custom-catalogue.so`` exists, that Arbor can load it, and
    that every mechanism in ``mechanisms/my_catalogue/*.mod`` is present.
    Rebuild after changing ``.mod`` files, upgrading Arbor, or changing
    OS/compiler; do not copy a ``.so`` between machines.

    Parameters
    ----------
    path : Path, optional
        Catalogue file. Default: ``CUSTOM_CATALOGUE_PATH``.

    Returns
    -------
    arbor.catalogue
    """
    catalogue_path = Path(path) if path is not None else CUSTOM_CATALOGUE_PATH
    if not catalogue_path.is_file():
        raise FileNotFoundError(
            f"Custom mechanism catalogue not found at {catalogue_path}. "
            "Build it from the repository root (see README Getting started):\n"
            f"  {_CATALOGUE_BUILD_CMD}"
        )
    logger.info("Loading custom catalogue from %s", catalogue_path)
    try:
        catalogue = A.load_catalogue(str(catalogue_path))
    except Exception as exc:
        raise RuntimeError(
            f"Could not load NMODL catalogue at {catalogue_path}. "
            "Rebuild after changing .mod files, upgrading Arbor, or changing "
            "OS/compiler. Do not copy a .so between machines.\n"
            f"  {_CATALOGUE_BUILD_CMD}"
        ) from exc
    expected = required_catalogue_mechanisms()
    missing = [
        name for name in expected if not _catalogue_has_mechanism(catalogue, name)
    ]
    if missing:
        raise RuntimeError(
            f"Catalogue at {catalogue_path} is missing mechanisms {missing}. "
            f"Expected {expected}. Rebuild with:\n"
            f"  {_CATALOGUE_BUILD_CMD}"
        )
    return catalogue

prepare_gap_junctions(swc_file, parameters)

Build gap junctions from CYCLE_BREAK and MULTI_NECK reconnect headers.

Returns:

Type	Description
dict[str, GapJunctionPoint]	Labels gj_0, gj_1, … in header order.

Examples:

>>> gjs = prepare_gap_junctions(swc_path, parameters)
>>> gjs["gj_0"].index_pair
(12, 13)

Source code in toric_spines_sim/model/gj.py

def prepare_gap_junctions(
    swc_file: Path,
    parameters: ParameterSet,
):
    """Build gap junctions from CYCLE_BREAK and MULTI_NECK reconnect headers.

    Returns
    -------
    dict[str, GapJunctionPoint]
        Labels ``gj_0``, ``gj_1``, … in header order.

    Examples
    --------
    >>> gjs = prepare_gap_junctions(swc_path, parameters)  # doctest: +SKIP
    >>> gjs["gj_0"].index_pair
    (12, 13)
    """
    weight = parameters["gj_weight"]
    reconnect_pairs = parse_reconnect_pairs(swc_file)
    gap_junctions = {}
    points_by_id = read_swc_points(swc_file)
    logger.debug(
        "Preparing gap junctions for %d reconnect pairs from %s, weight=%f",
        len(reconnect_pairs),
        swc_file,
        weight,
    )
    for n, (i, j) in enumerate(reconnect_pairs):
        if i not in points_by_id or j not in points_by_id:
            logger.error("Node %s or %s not found in SWC file %s", i, j, swc_file)
            raise ValueError(f"Node {i} or {j} not found in SWC file")
        xi, yi, zi, _ri = points_by_id[i]
        xj, yj, zj, _rj = points_by_id[j]
        loc_a = (xi, yi, zi)
        loc_b = (xj, yj, zj)
        gap_junctions[f"gj_{n}"] = GapJunctionPoint(
            (i, j),
            loc_a,
            weight,
            location_b=None if loc_b == loc_a else loc_b,
        )
    logger.debug(
        "Prepared %d gap junctions from %s weight=%f",
        len(gap_junctions),
        swc_file,
        weight,
    )
    return gap_junctions