[WIP] 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]
+    v : 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.
+    v : float
+        start 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, v=-65, end_pts=None):
         self._L = L
         self._diam = diam
         self._Ra = Ra
         self._cm = cm
+        self._v = v  # initial voltage
         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}, v={self.v}"
 
     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["v"] = self._v
         # 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 v(self):
+        return self._v
+
     @property
     def end_pts(self):
         return self._end_pts
@@ -574,6 +584,7 @@ def _set_biophysics(self, sections):
             sec = self._nrn_sections[sec_name]
             for mech_name, p_mech in section.mechs.items():
                 sec.insert(mech_name)
+                setattr(sec, "v", section.v)
                 for attr, val in p_mech.items():
                     if isinstance(val, list):
                         seg_xs, seg_vals = val[0], val[1]
@@ -642,6 +653,7 @@ 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
+            sec.v = sections[sec_name].v
 
         if cell_tree is None:
             cell_tree = dict()
@@ -1062,7 +1074,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, v=None):
         """Change attributes of section specified by `sec_name`
 
         Parameters
@@ -1077,6 +1089,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.
+        v : float | int | None
+            start value for membrane potential in millivolts.
 
         Notes
         -----
@@ -1101,4 +1115,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 v is not None:
+            _validate_type(v, (float, int), "v")
+            self.sections[sec_name]._v = v
+
         self._update_end_pts()

hnn_core/cells_default.py

@@ -15,7 +15,9 @@
 # units for taur: ms
 
 
-def _get_dends(params, cell_type, section_names):
+# KD: comment: initialize membrane potential here as it's not overridden by
+# h.finitialize unless called as h.finitialize(-65)
+def _get_dends(params, cell_type, section_names, v_init={"all": -65}):
     """Convert a flat dictionary to a nested dictionary.
 
     Returns
@@ -34,22 +36,64 @@ def _get_dends(params, cell_type, section_names):
                 # map apicaltrunk -> apical_trunk etc.
                 middle = section_name.replace("_", "")
             dend_prop[key] = params[f"{cell_type}_{middle}_{key}"]
+            if len(v_init) == 1:
+                v = v_init["all"]
+            else:
+                v = v_init[section_name]
+
         sections[section_name] = Section(
             L=dend_prop["L"],
             diam=dend_prop["diam"],
             Ra=dend_prop["Ra"],
             cm=dend_prop["cm"],
+            v=v,
         )
     return sections
 
 
-def _get_pyr_soma(p_all, cell_type):
+# KD: In the new model, the basal dendrites are differently tuned from the apical
+# dendrites.
+def _get_basal(params, cell_type, section_names, v_init={"all": -65}):
+    """Convert a flat dictionary to a nested dictionary.
+
+    Returns
+    -------
+    sections : dict
+        Dictionary of sections. Keys are section names
+    """
+    prop_names = ["L", "diam", "Ra", "cm"]
+    sections = dict()
+    for section_name in section_names:
+        dend_prop = dict()
+        middle = section_name.replace("_", "")
+        for key in prop_names:
+            if key in ["Ra", "cm"]:
+                middle = "basal"
+            else:
+                # map apicaltrunk -> apical_trunk etc.
+                middle = section_name.replace("_", "")
+            dend_prop[key] = params[f"{cell_type}_{middle}_{key}"]
+            if len(v_init) == 1:
+                v = v_init["all"]
+            else:
+                v = v_init[section_name]
+        sections[section_name] = Section(
+            L=dend_prop["L"],
+            diam=dend_prop["diam"],
+            Ra=dend_prop["Ra"],
+            cm=dend_prop["cm"],
+            v=v,
+        )
+
+
+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"],
+        v=v_init,
     )
 
 
@@ -60,6 +104,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 +116,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 +210,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 +308,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,
+        v=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 "v" in section_data.keys():
+                sections[section_name] = Section(
+                    L=section_data["L"],
+                    diam=section_data["diam"],
+                    cm=section_data["cm"],
+                    Ra=section_data["Ra"],
+                    v=section_data["v"],
+                    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.