[MRG] feat: Extracting Katharina's voltage-initialization changes

hnn_core/cell.py

@@ -180,7 +180,9 @@ class Section:
         membrane capacitance in micro-Farads.
     Ra : float
         axial resistivity in ohm-cm
-    end_pts : list of [x, y, z]
+    v0 : float
+        start value for membrane potential in millivolts.
+    end_pts : list of [x, y, z], optional
         The start and stop points of the section.
 
     Attributes
@@ -201,15 +203,18 @@ class Section:
         membrane capacitance in micro-Farads.
     Ra : float
         axial resistivity in ohm-cm.
+    v0 : float
+        Initial value for membrane potential in millivolts.
     nseg : int
         Number of segments in the section
     """
 
-    def __init__(self, L, diam, Ra, cm, end_pts=None):
+    def __init__(self, L, diam, Ra, cm, v0=-65, end_pts=None):
         self._L = L
         self._diam = diam
         self._Ra = Ra
         self._cm = cm
+        self._v0 = v0
         if end_pts is None:
             end_pts = list()
         self._end_pts = end_pts
@@ -221,7 +226,7 @@ def __init__(self, L, diam, Ra, cm, end_pts=None):
         self.nseg = _get_nseg(self.L)
 
     def __repr__(self):
-        return f"L={self.L}, diam={self.diam}, cm={self.cm}, Ra={self.Ra}"
+        return f"L={self.L}, diam={self.diam}, cm={self.cm}, Ra={self.Ra}, v0={self.v0}"
 
     def __eq__(self, other):
         if not isinstance(other, Section):
@@ -266,6 +271,7 @@ def to_dict(self):
         section_data["Ra"] = self.Ra
         section_data["end_pts"] = self.end_pts
         section_data["nseg"] = self.nseg
+        section_data["v0"] = self.v0
         # Need to solve the partial function problem
         # in mechs
         section_data["mechs"] = self.mechs
@@ -288,6 +294,10 @@ def cm(self):
     def Ra(self):
         return self._Ra
 
+    @property
+    def v0(self):
+        return self._v0
+
     @property
     def end_pts(self):
         return self._end_pts
@@ -574,6 +584,9 @@ def _set_biophysics(self, sections):
             sec = self._nrn_sections[sec_name]
             for mech_name, p_mech in section.mechs.items():
                 sec.insert(mech_name)
+                # This is one of the two places where we are actually applying our
+                # initial voltage
+                setattr(sec, "v", section.v0)
                 for attr, val in p_mech.items():
                     if isinstance(val, list):
                         seg_xs, seg_vals = val[0], val[1]
@@ -642,6 +655,9 @@ def _create_sections(self, sections, cell_tree):
             sec.Ra = sections[sec_name].Ra
             sec.cm = sections[sec_name].cm
             sec.nseg = sections[sec_name].nseg
+            # This is one of the two places where we are actually applying our initial
+            # voltage
+            sec.v = sections[sec_name].v0
 
         if cell_tree is None:
             cell_tree = dict()
@@ -1062,7 +1078,7 @@ def _update_end_pts(self):
         # of sections.
         self.define_shape(("soma", 0))
 
-    def modify_section(self, sec_name, L=None, diam=None, cm=None, Ra=None):
+    def modify_section(self, sec_name, L=None, diam=None, cm=None, Ra=None, v0=None):
         """Change attributes of section specified by `sec_name`
 
         Parameters
@@ -1077,6 +1093,8 @@ def modify_section(self, sec_name, L=None, diam=None, cm=None, Ra=None):
             membrane capacitance in micro-Farads.
         Ra : float | int | None
             axial resistivity in ohm-cm.
+        v0 : float | int | None
+            start value for membrane potential in millivolts.
 
         Notes
         -----
@@ -1101,4 +1119,8 @@ def modify_section(self, sec_name, L=None, diam=None, cm=None, Ra=None):
             _validate_type(Ra, (float, int), "Ra")
             self.sections[sec_name]._Ra = Ra
 
+        if v0 is not None:
+            _validate_type(v0, (float, int), "v0")
+            self.sections[sec_name]._v0 = v0
+
         self._update_end_pts()

hnn_core/cells_default.py

@@ -15,41 +15,101 @@
 # units for taur: ms
 
 
-def _get_dends(params, cell_type, section_names):
-    """Convert a flat dictionary to a nested dictionary.
+def _get_dends(
+        params,
+        cell_type,
+        section_names,
+        v_init={"all": -65},
+        is_basal_specific=False,
+):
+    """Create dendritic Section objects from flat parameter dictionary.
+
+    Extracts geometric and electrical properties (length, diameter, axial resistance,
+    membrane capacitance) from a flat parameter dictionary, takes initial membrane
+    voltage from its, and constructs Section objects for each dendritic
+    compartment. Handles parameter key name transformations (e.g., 'apical_trunk' ->
+    'apicaltrunk') required for lookup in the parameter dictionary.
+
+    Parameters
+    ----------
+    params : dict
+        Flat dictionary containing cell parameters with keys formatted as
+        '{cell_type}_{section}_{property}' (e.g., 'L5Pyr_apicaltrunk_L').
+        'Ra' and 'cm' use "dend" as the middle component rather than specific
+        section names. This 'params' dictionary is expected to be constructed using
+        functions like `params_default.py::get_L2Pyr_params_default`.
+    cell_type : str, {'L2Pyr', 'L5Pyr'}
+        Cell type identifier used as prefix in parameter key lookups.
+    section_names : list of str
+        Names of dendritic sections to create (e.g., ['apical_trunk',
+        'apical_1', 'basal_2']). Underscores are removed for parameter
+        lookups except for 'Ra' and 'cm'.
+    v_init : dict, default={"all": -65}
+        Initial membrane potential in mV. If dict contains single key "all",
+        that value is applied to all sections. Otherwise, keys must match
+        'section_names' for section-specific initialization.
+    is_basal_specific : bool, default=False
+        Flag indicating whether or not to use the (Duecker 2025) model's custom basal
+        dendrite parameters. If True, this will read the 'Ra' and 'cm' parameters from
+        'params' using '{cell_type}_basal_{property}' instead of the default
+        '{cell_type}_dend_{property}' naming scheme.
 
     Returns
     -------
     sections : dict
-        Dictionary of sections. Keys are section names
+        Dictionary mapping section names (str) to Section objects with attributes L,
+        diam, Ra, and cm set from 'params', and v0 set from argument.
+
+    Notes
+    -----
+    - KD: This function is where the initial voltages for the dendritic sections are
+      set; these voltages are not overridden by `h.finitialize` unless called with a
+      value, e.g. `h.finitialize(-65)`.
+    - The 'v0' (initial voltage) parameter is handled separately from other properties as
+      it is a newer addition not found in legacy parameter files.
+    - In the (Jones et al., 2009) model, this is used to construct both apical and basal
+      dendrite sections. In the newer (Duecker 2025) model, this is only used for the
+      apical dendrite sections.
     """
     prop_names = ["L", "diam", "Ra", "cm"]
     sections = dict()
     for section_name in section_names:
         dend_prop = dict()
         for key in prop_names:
             if key in ["Ra", "cm"]:
-                middle = "dend"
+                if is_basal_specific:
+                    middle = "basal"
+                else:
+                    middle = "dend"
             else:
                 # map apicaltrunk -> apical_trunk etc.
                 middle = section_name.replace("_", "")
             dend_prop[key] = params[f"{cell_type}_{middle}_{key}"]
+        # v0 is handled separately since it is "newer", and will never be found in the
+        # `params` input.
+        if len(v_init) == 1:
+            dend_prop["v0"] = v_init["all"]
+        else:
+            dend_prop["v0"] = v_init[section_name]
+
         sections[section_name] = Section(
             L=dend_prop["L"],
             diam=dend_prop["diam"],
             Ra=dend_prop["Ra"],
             cm=dend_prop["cm"],
+            v0=dend_prop["v0"],
         )
     return sections
 
 
-def _get_pyr_soma(p_all, cell_type):
+def _get_pyr_soma(p_all, cell_type, v_init=-65):
     """Get somatic properties."""
     return Section(
         L=p_all[f"{cell_type}_soma_L"],
         diam=p_all[f"{cell_type}_soma_diam"],
         cm=p_all[f"{cell_type}_soma_cm"],
         Ra=p_all[f"{cell_type}_soma_Ra"],
+        v0=v_init,
     )
 
 
@@ -60,6 +120,8 @@ def _cell_L2Pyr(override_params, pos=(0.0, 0.0, 0), gid=0.0):
         assert isinstance(override_params, dict)
         p_all = compare_dictionaries(p_all, override_params)
 
+    all_v_init = -71.46
+
     section_names = [
         "apical_trunk",
         "apical_1",
@@ -70,8 +132,19 @@ def _cell_L2Pyr(override_params, pos=(0.0, 0.0, 0), gid=0.0):
         "basal_3",
     ]
 
-    sections = _get_dends(p_all, cell_type="L2Pyr", section_names=section_names)
-    sections["soma"] = _get_pyr_soma(p_all, "L2Pyr")
+    sections = _get_dends(
+        p_all,
+        cell_type="L2Pyr",
+        section_names=section_names,
+        v_init={
+            "all": all_v_init,
+        },
+    )
+    sections["soma"] = _get_pyr_soma(
+        p_all,
+        "L2Pyr",
+        v_init=all_v_init,
+    )
 
     end_pts = {
         "soma": [[-50, 0, 765], [-50, 0, 778]],
@@ -153,8 +226,29 @@ def _cell_L5Pyr(override_params, pos=(0.0, 0.0, 0), gid=0.0):
         "basal_3",
     ]
 
-    sections = _get_dends(p_all, cell_type="L5Pyr", section_names=section_names)
-    sections["soma"] = _get_pyr_soma(p_all, "L5Pyr")
+    v_init = {
+        "apical_1": -71.32,
+        "apical_2": -69.08,
+        "apical_tuft": -67.30,
+        "apical_trunk": -72,
+        "soma": -72.0,
+        "basal_1": -72,
+        "basal_2": -72,
+        "basal_3": -72,
+        "apical_oblique": -72,
+    }
+
+    sections = _get_dends(
+        p_all,
+        cell_type="L5Pyr",
+        section_names=section_names,
+        v_init=v_init,
+    )
+    sections["soma"] = _get_pyr_soma(
+        p_all,
+        "L5Pyr",
+        v_init=-72,
+    )
 
     end_pts = {
         "soma": [[0, 0, 0], [0, 0, 23]],
@@ -230,9 +324,16 @@ def _cell_L5Pyr(override_params, pos=(0.0, 0.0, 0), gid=0.0):
     )
 
 
-def _get_basket_soma(cell_name):
+def _get_basket_soma(cell_name, v_init=-64.9737):
     end_pts = [[0, 0, 0], [0, 0, 39.0]]
-    return Section(L=39.0, diam=20.0, cm=0.85, Ra=200.0, end_pts=end_pts)
+    return Section(
+        L=39.0,
+        diam=20.0,
+        cm=0.85,
+        Ra=200.0,
+        v0=v_init,
+        end_pts=end_pts,
+    )
 
 
 def _get_pyr_syn_props(p_all, cell_type):

hnn_core/hnn_io.py

@@ -166,13 +166,25 @@ def _read_cell_types(cell_types_data):
         sections_data = cell_data["sections"]
         for section_name in sections_data:
             section_data = sections_data[section_name]
-            sections[section_name] = Section(
-                L=section_data["L"],
-                diam=section_data["diam"],
-                cm=section_data["cm"],
-                Ra=section_data["Ra"],
-                end_pts=section_data["end_pts"],
-            )
+            if "v0" in section_data.keys():
+                sections[section_name] = Section(
+                    L=section_data["L"],
+                    diam=section_data["diam"],
+                    cm=section_data["cm"],
+                    Ra=section_data["Ra"],
+                    v0=section_data["v0"],
+                    end_pts=section_data["end_pts"],
+                )
+            else:
+                # Yet more legacy backwards-compatibility
+                sections[section_name] = Section(
+                    L=section_data["L"],
+                    diam=section_data["diam"],
+                    cm=section_data["cm"],
+                    Ra=section_data["Ra"],
+                    end_pts=section_data["end_pts"],
+                )
+
             # Set section attributes
             sections[section_name].syns = section_data["syns"]
             sections[section_name].mechs = section_data["mechs"]

hnn_core/network_builder.py

@@ -371,8 +371,6 @@ def _build(self):
             record_ca=record_ca,
         )
 
-        self.state_init()
-
         # set to record spikes, somatic voltages, and extracellular potentials
         self._spike_times = h.Vector()
         self._spike_gids = h.Vector()
@@ -649,37 +647,6 @@ def aggregate_data(self, n_samples):
 
         _PC.barrier()  # get all nodes to this place before continuing
 
-    def state_init(self):
-        """Initializes the state closer to baseline."""
-
-        for cell in self._cells:
-            seclist = h.SectionList()
-            seclist.wholetree(sec=cell._nrn_sections["soma"])
-            src_type = self.net.gid_to_type(cell.gid)
-            cell_metadata = self.net.cell_types[src_type]["cell_metadata"]
-            # initializing segment voltages from cell_metadata
-            for sect in seclist:
-                for seg in sect:
-                    if (
-                        cell_metadata.get("morpho_type") == "pyramidal"
-                        and cell_metadata.get("layer") == "2"
-                    ):
-                        seg.v = -71.46
-                    elif (
-                        cell_metadata.get("morpho_type") == "pyramidal"
-                        and cell_metadata.get("layer") == "5"
-                    ):
-                        if sect.name() == f"{_short_name(src_type)}_apical_1":
-                            seg.v = -71.32
-                        elif sect.name() == f"{_short_name(src_type)}_apical_2":
-                            seg.v = -69.08
-                        elif sect.name() == f"{_short_name(src_type)}_apical_tuft":
-                            seg.v = -67.30
-                        else:
-                            seg.v = -72.0
-                    elif cell_metadata.get("morpho_type") == "basket":
-                        seg.v = -64.9737
-
     def _clear_neuron_objects(self):
         """Clear up NEURON internal gid and reference information.