[MRG] add method for updating cell biophysics from Network instance

doc/Makefile

@@ -2,7 +2,7 @@
 #
 
 # You can set these variables from the command line.
-SPHINXOPTS    =
+SPHINXOPTS    = -v
 SPHINXBUILD   = sphinx-build
 PAPER         =
 BUILDDIR      = _build

doc/whats_new.rst

@@ -20,6 +20,8 @@ Changelog
 
 - Compute dipole component in z-direction automatically from cell morphology instead of hard coding, by `Mainak Jas`_ in  `#327 <https://github.com/jonescompneurolab/hnn-core/pull/320>`_
 
+- Store :class:`~hnn_core.Cell` instances in :class:`~hnn_core.Network`'s :attr:`~/hnn_core.Network.cells` attribute by `Ryan Thorpe`_ in `#321 <https://github.com/jonescompneurolab/hnn-core/pull/321>`_
+
 Bug
 ~~~
 
@@ -30,7 +32,10 @@ API
 - Remove :class:`~hnn_core.L2Pyr`, :class:`~hnn_core.L5Pyr`, :class:`~hnn_core.L2Basket`, and :class:`~hnn_core.L5Basket` classes
   in favor of instantation through functions and a more consistent :class:`~hnn_core.Cell` class by `Mainak Jas`_ in  `#322 <https://github.com/jonescompneurolab/hnn-core/pull/320>`_
 
-- Remove parameter `distribution` in :func:`~hnn_core.Network.add_bursty_drive`. The distribution is now Gaussian by default, by `Mainak Jas` in `#330 <https://github.com/jonescompneurolab/hnn-core/pull/330>`_
+- Remove parameter `distribution` in :func:`~hnn_core.Network.add_bursty_drive`. The distribution is now Gaussian by default, by `Mainak Jas`_ in `#330 <https://github.com/jonescompneurolab/hnn-core/pull/330>`_
+
+- New API for accessing and modifying :class:`~hnn_core.Cell` attributes (e.g., synapse and biophysics parameters) as cells are now instantiated from template cells specified
+  in a :class:`~hnn_core.Network` instance's :attr:`~/hnn_core.Network.cell_types` attribute by `Ryan Thorpe`_ in `#321 <https://github.com/jonescompneurolab/hnn-core/pull/321>`_
 
 .. _0.1:
 
@@ -74,9 +79,9 @@ Changelog
 
 - Add functions for plotting power spectral density (:func:`~hnn_core.viz.plot_psd`) and Morlet time-frequency representations (:func:`~hnn_core.viz.plot_tfr_morlet`), by `Christopher Bailey`_ in `#264 <https://github.com/jonescompneurolab/hnn-core/pull/264>`_
 
-- Add y-label units (nAm) to all visualisation functions involving dipole moments, by `Christopher Bailey`_ in `#264 <https://github.com/jonescompneurolab/hnn-core/pull/264>`_ 
+- Add y-label units (nAm) to all visualisation functions involving dipole moments, by `Christopher Bailey`_ in `#264 <https://github.com/jonescompneurolab/hnn-core/pull/264>`_
 
-- Add Savitzky-Golay filtering method :meth:`~hnn_core.dipole.Dipole.savgol_filter` to ``Dipole``; copied from ``mne-python`` :meth:`~mne.Evoked.savgol_filter`, by `Christopher Bailey`_ in `#264 <https://github.com/jonescompneurolab/hnn-core/pull/264>`_ 
+- Add Savitzky-Golay filtering method :meth:`~hnn_core.dipole.Dipole.savgol_filter` to ``Dipole``; copied from ``mne-python`` :meth:`~mne.Evoked.savgol_filter`, by `Christopher Bailey`_ in `#264 <https://github.com/jonescompneurolab/hnn-core/pull/264>`_
 
 Bug
 ~~~

examples/plot_simulate_gamma.py

@@ -82,7 +82,7 @@
 dpls[trial_idx].plot_tfr_morlet(freqs, n_cycles=7, tmin=tmin, ax=axes[1])
 
 ###############################################################################
-# As a final exercise, let us try to re-run the simulation with a tonic bias
+# Now, let us try to re-run the simulation with a tonic bias
 # applied to the L5 Pyramidal cells. Notice that the oscillation waveform is
 # more regular, with less noise due to the fact that the tonic depolarization
 # dominates over the influence of the Poisson drive. By default, a tonic bias
@@ -108,6 +108,21 @@
 from hnn_core.viz import plot_psd
 plot_psd(dpls[trial_idx], fmin=20., fmax=100., tmin=tmin)
 
+###############################################################################
+# Finally, we demonstrate the mechanistic link between PING and the GABAA decay
+# time constant (`tau2`). Using the same network/drive configuration as before,
+# we decrease `tau2` from 5 to 2 ms. This will shorten the effective
+# refactory period between L5 pyramidal cell spikes and increase the PING
+# frequency from ~50 to ~65 Hz.
+net.cell_types['L5_pyramidal'].p_syn['gabaa']['tau2'] = 2
+dpls = simulate_dipole(net, n_trials=1)
+
+fig, axes = plt.subplots(3, 1, sharex=True, figsize=(6, 6),
+                         constrained_layout=True)
+dpls[trial_idx].plot(ax=axes[0], show=False)
+net.cell_response.plot_spikes_raster(ax=axes[1], show=False)
+dpls[trial_idx].plot_tfr_morlet(freqs, n_cycles=7, tmin=tmin, ax=axes[2])
+
 ###############################################################################
 # References
 # ----------

hnn_core/cell.py

@@ -91,6 +91,29 @@ class Cell:
         The name of the cell.
     pos : tuple
         The (x, y, z) coordinates.
+    p_secs : nested dict
+        Dictionary with keys as section name.
+        p_secs[sec_name] is a dictionary with keys
+        L, diam, Ra, cm, syns and mech.
+        syns is a list specifying the synapses at that section.
+        The properties of syn are specified in p_syn.
+        mech is a dict with keys as the mechanism names. The
+        values are dictionaries with properties of the mechanism.
+    p_syn : dict of dict
+        Keys are name of synaptic mechanism. Each synaptic mechanism
+        has keys for parameters of the mechanism, e.g., 'e', 'tau1',
+        'tau2'.
+    topology : list of list
+        The topology of cell sections. Each element is a list of
+        4 items in the format
+        [parent_sec, parent_loc, child_sec, child_loc] where
+        parent_sec and parent_loc are float between 0 and 1
+        specifying the location in the section to connect and
+        parent_sec and child_sec are names of the connecting
+        sections.
+    sect_loc : dict of list
+        Can have keys 'proximal' or 'distal' each containing
+        names of section locations that are proximal or distal.
     gid : int or None (optional)
         Each cell in a network is uniquely identified by it's "global ID": GID.
         The GID is an integer from 0 to n_cells, or None if the cell is not
@@ -122,15 +145,37 @@ class Cell:
     sect_loc : dict of list
         Can have keys 'proximal' or 'distal' each containing
         names of section locations that are proximal or distal.
+
+    Examples
+    --------
+    >>> p_secs = {
+        'soma':
+        {
+            'L': 39,
+            'diam': 20,
+            'cm': 0.85,
+            'Ra': 200.,
+            'sec_pts': [[0, 0, 0], [0, 39., 0]],
+            'syns': ['ampa', 'gabaa', 'nmda'],
+            'mechs' : {
+                'ca': {
+                    'gbar_ca': 60
+                }
+            }
+        }
+    }
     """
 
-    def __init__(self, name, pos, gid=None):
+    def __init__(self, name, pos, p_secs, p_syn, topology, sect_loc, gid=None):
         self.name = name
         self.pos = pos
+        self.p_secs = p_secs
+        self.p_syn = p_syn
+        self.topology = topology
+        self.sect_loc = sect_loc
         self.sections = dict()
         self.synapses = dict()
         self.dipole_pp = list()
-        self.sect_loc = dict()
         self.rec_v = h.Vector()
         self.rec_i = dict()
         # insert iclamp
@@ -232,60 +277,27 @@ def _create_sections(self, p_secs, topology):
             child_loc = connection[3]
             child_sec.connect(parent_sec, parent_loc, child_loc)
 
-    def build(self, p_secs, p_syn, topology, sect_loc):
-        """Build cell in Neuron.
+    def build(self, sec_name_apical=None):
+        """Build cell in Neuron and insert dipole if applicable.
 
         Parameters
         ----------
-        p_secs : nested dict
-            Dictionary with keys as section name.
-            p_secs[sec_name] is a dictionary with keys
-            L, diam, Ra, cm, syns and mech.
-            syns is a list specifying the synapses at that section.
-            The properties of syn are specified in p_syn.
-            mech is a dict with keys as the mechanism names. The
-            values are dictionaries with properties of the mechanism.
-        p_syn : dict of dict
-            Keys are name of synaptic mechanism. Each synaptic mechanism
-            has keys for parameters of the mechanism, e.g., 'e', 'tau1',
-            'tau2'.
-        topology : list of list
-            The topology of cell sections. Each element is a list of
-            4 items in the format
-            [parent_sec, parent_loc, child_sec, child_loc] where
-            parent_sec and parent_loc are float between 0 and 1
-            specifying the location in the section to connect and
-            parent_sec and child_sec are names of the connecting
-            sections.
-        sect_loc : dict of list
-            Can have keys 'proximal' or 'distal' each containing
-            names of section locations that are proximal or distal.
-
-        Examples
-        --------
-        >>> p_secs = {
-            'soma':
-            {
-                'L': 39,
-                'diam': 20,
-                'cm': 0.85,
-                'Ra': 200.,
-                'sec_pts': [[0, 0, 0], [0, 39., 0]],
-                'syns': ['ampa', 'gabaa', 'nmda'],
-                'mechs' : {
-                    'ca': {
-                        'gbar_ca': 60
-                    }
-                }
-            }
-        }
+        sec_name_apical : str | None
+            If not None, a dipole will be inserted in this cell in alignment
+            with this section. The section should belong to the apical dendrite
+            of a pyramidal neuron.
         """
-        if 'soma' not in p_secs:
+        if 'soma' not in self.p_secs:
             raise KeyError('soma must be defined for cell')
-        self._create_sections(p_secs, topology)
-        self._create_synapses(p_secs, p_syn)
-        self._set_biophysics(p_secs)
-        self.sect_loc = sect_loc
+        self._create_sections(self.p_secs, self.topology)
+        self._create_synapses(self.p_secs, self.p_syn)
+        self._set_biophysics(self.p_secs)
+        if sec_name_apical in self.sections:
+            self._insert_dipole(sec_name_apical)
+        elif sec_name_apical is not None:
+            raise ValueError(f'sec_name_apical must be an existing '
+                             f'section of the current cell or None. '
+                             f'Got {sec_name_apical}.')
 
     def move_to_pos(self):
         """Move cell to position."""
@@ -306,28 +318,20 @@ def move_to_pos(self):
     # 2. a list needs to be created with a Dipole (Point Process) in each
     #    section at position 1
     # In Cell() and not Pyr() for future possibilities
-    def insert_dipole(self, p_secs, sec_name_apical):
+    def _insert_dipole(self, sec_name_apical):
         """Insert dipole into each section of this cell.
 
         Parameters
         ----------
-        p_secs : nested dict
-            Dictionary with keys as section name.
-            p_secs[sec_name] is a dictionary with keys
-            L, diam, Ra, cm, syns and mech.
-            syns is a list specifying the synapses at that section.
-            The properties of syn are specified in p_syn.
-            mech is a dict with keys as the mechanism names. The
-            values are dictionaries with properties of the mechanism.
         sec_name_apical : str
             The name of the section along which dipole moment is calculated.
         """
         self.dpl_vec = h.Vector(1)
         self.dpl_ref = self.dpl_vec._ref_x[0]
-        cos_thetas = _get_cos_theta(p_secs, 'apical_trunk')
+        cos_thetas = _get_cos_theta(self.p_secs, 'apical_trunk')
 
         # setting pointers and ztan values
-        for sect_name in p_secs:
+        for sect_name in self.p_secs:
             sect = self.sections[sect_name]
             sect.insert('dipole')
 

hnn_core/cell_response.py

@@ -27,6 +27,12 @@ class CellResponse(object):
         The inner list contains the type of spike (e.g., evprox1
         or L2_pyramidal) that occured at the corresonding time stamp.
         Each gid corresponds to a type via Network().gid_ranges.
+    times : numpy array | None
+        Array of time points for samples in continuous data.
+        This includes vsoma and isoma.
+    cell_type_names : list
+        List of unique cell type names that are explicitly modeled in the
+        network
 
     Attributes
     ----------
@@ -48,7 +54,6 @@ class CellResponse(object):
     isoma : list (n_trials,) of dict, shape
         Each element of the outer list is a trial.
         Dictionary indexed by gids containing somatic currents.
-
     times : numpy array
         Array of time points for samples in continuous data.
         This includes vsoma and isoma.
@@ -70,7 +75,10 @@ class CellResponse(object):
     """
 
     def __init__(self, spike_times=None, spike_gids=None, spike_types=None,
-                 times=None):
+                 times=None, cell_type_names=['L2_basket',
+                                              'L2_pyramidal',
+                                              'L5_basket',
+                                              'L5_pyramidal']):
         if spike_times is None:
             spike_times = list()
         if spike_gids is None:
@@ -106,6 +114,7 @@ def __init__(self, spike_times=None, spike_gids=None, spike_types=None,
             if not isinstance(times, np.ndarray):
                 raise TypeError("'times' is an np.ndarray of simulation times")
         self._times = times
+        self._cell_type_names = cell_type_names
 
     def __repr__(self):
         class_name = self.__class__.__name__
@@ -281,7 +290,6 @@ def mean_rates(self, tstart, tstop, gid_ranges, mean_type='all'):
         spike_rate : dict
             Dictionary with keys 'L5_pyramidal', 'L5_basket', etc.
         """
-        cell_types = ['L5_pyramidal', 'L5_basket', 'L2_pyramidal', 'L2_basket']
         spike_rates = dict()
 
         if mean_type not in ['all', 'trial', 'cell']:
@@ -295,7 +303,7 @@ def mean_rates(self, tstart, tstop, gid_ranges, mean_type='all'):
         elif tstop <= tstart:
             raise ValueError('tstop must be greater than tstart')
 
-        for cell_type in cell_types:
+        for cell_type in self._cell_type_names:
             cell_type_gids = np.array(gid_ranges[cell_type])
             n_trials, n_cells = len(self._spike_times), len(cell_type_gids)
             gid_spike_rate = np.zeros((n_trials, n_cells))

hnn_core/cells_default.py

@@ -146,15 +146,20 @@ def _get_mechanisms(p_all, cell_type, section_names, mechanisms):
     return mech_props
 
 
-def basket(pos, cell_name='L2Basket', gid=None):
+def _set_variable_mech(dist_from_soma):
+    """Set a cell mechanism based on its distance from the soma"""
+    return 1e-6 * np.exp(3e-3 * dist_from_soma)
+
+
+def basket(cell_name, pos=(0, 0, 0), gid=None):
     """Get layer 2 / layer 5 basket cells.
 
     Parameters
     ----------
-    pos : tuple
-        Coordinates of cell soma in xyz-space
     cell_name : str
         The name of the cell.
+    pos : tuple
+        Coordinates of cell soma in xyz-space
     gid : int or None (optional)
         Each cell in a network is uniquely identified by it's "global ID": GID.
         The GID is an integer from 0 to n_cells, or None if the cell is not
@@ -182,21 +187,23 @@ def basket(pos, cell_name='L2Basket', gid=None):
     for sec_name in p_secs:
         p_secs[sec_name]['sec_pts'] = sec_pts[sec_name]
 
-    cell = Cell(cell_name, pos=pos, gid=gid)
-    cell.build(p_secs, p_syn, topology, sect_loc)
-
-    return cell
+    return Cell(cell_name, pos,
+                p_secs=p_secs,
+                p_syn=p_syn,
+                topology=topology,
+                sect_loc=sect_loc,
+                gid=gid)
 
 
-def pyramidal(pos, cell_name, override_params=None, gid=None):
+def pyramidal(cell_name, pos=(0, 0, 0), override_params=None, gid=None):
     """Pyramidal neuron.
 
     Parameters
     ----------
-    pos : tuple
-        Coordinates of cell soma in xyz-space
     cell_name : str
         'L5Pyr' or 'L2Pyr'. The pyramidal cell type.
+    pos : tuple
+        Coordinates of cell soma in xyz-space
     override_params : dict or None (optional)
         Parameters specific to L2 pyramidal neurons to override the default set
     gid : int or None (optional)
@@ -258,7 +265,7 @@ def pyramidal(pos, cell_name, override_params=None, gid=None):
             syns = list(p_syn.keys())
             if cell_name == 'L5Pyr':
                 p_secs[key]['mechs'][
-                    'ar']['gbar_ar'] = lambda x: 1e-6 * np.exp(3e-3 * x)
+                    'ar']['gbar_ar'] = _set_variable_mech
         p_secs[key]['syns'] = syns
 
     for sec_name in p_secs:
@@ -267,10 +274,9 @@ def pyramidal(pos, cell_name, override_params=None, gid=None):
     sect_loc = {'proximal': ['apical_oblique', 'basal_2', 'basal_3'],
                 'distal': ['apical_tuft']}
 
-    cell = Cell(name=cell_name, pos=pos, gid=gid)
-    cell.build(p_secs, p_syn, topology, sect_loc=sect_loc)
-
-    # insert dipole
-    cell.insert_dipole(p_secs, 'apical_trunk')
-
-    return cell
+    return Cell(cell_name, pos,
+                p_secs=p_secs,
+                p_syn=p_syn,
+                topology=topology,
+                sect_loc=sect_loc,
+                gid=gid)