[MRG] Add functionality to record somatic voltage from all cells in simulation

doc/whats_new.rst

@@ -31,6 +31,8 @@ Changelog
 
 - Add function to compute mean spike rates with user specified calculation type, by `Nick Tolley`_ and `Mainak Jas`_ in `#155 <https://github.com/jonescompneurolab/hnn-core/pull/155>`_
 
+- Add ability to record somatic voltages from all cells, by `Nick Tolley`_ in `#190 <https://github.com/jonescompneurolab/hnn-core/pull/190>`_
+
 Bug
 ~~~
 

examples/plot_firing_pattern.py

@@ -15,8 +15,7 @@
 # Let us import hnn_core
 
 import hnn_core
-from hnn_core import read_params, Network
-from hnn_core.network_builder import NetworkBuilder
+from hnn_core import read_params, Network, simulate_dipole
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 
@@ -26,30 +25,53 @@
 params = read_params(params_fname)
 
 ###############################################################################
-# Now let's build the network
+# Now let's build the network with somatic voltage recordings enabled
 import matplotlib.pyplot as plt
 
 net = Network(params)
-with NetworkBuilder(net) as network_builder:
-    network_builder.cells[0].plot_voltage()
+dpls = simulate_dipole(net, n_trials=1, record_vsoma=True)
 
-    # The cells are stored in the network object as a list
-    cells = network_builder.cells
-    print(cells[:5])
+###############################################################################
+# The cell IDs (gids) are stored in the network object as a dictionary
+gid_dict = net.gid_dict
+print(net.gid_dict)
 
-    # We have different kinds of cells with different cell IDs (gids)
-    gids = [0, 35, 135, 170]
-    for gid in gids:
-        print(cells[gid].name)
+###############################################################################
+# Simulated cell responses are stored in the Spikes object as a dictionary.
+trial_idx = 0
+vsoma = net.spikes.vsoma[trial_idx]
+print(vsoma.keys())
 
-    # We can plot the firing pattern of individual cells
-    network_builder.cells[0].plot_voltage()
-    plt.title('%s (gid=%d)' % (cells[0].name, gid))
+###############################################################################
+# We can plot the firing pattern of individual cells by indexing with the gid
+gid = 170
+times = net.spikes.times[trial_idx]
+plt.figure(figsize=(4, 4))
+plt.plot(times, vsoma[gid])
+plt.title('%s (gid=%d)' % (net.gid_to_type(gid), gid))
+plt.xlabel('Time (ms)')
+plt.ylabel('Voltage (mV)')
+plt.show()
 
 ###############################################################################
-# Let's do this for the rest of the cell types with a new NetworkBuilder object
-with NetworkBuilder(net) as network_builder:
-    fig, axes = plt.subplots(1, 2, sharey=True, figsize=(8, 4))
-    for gid, ax in zip([35, 170], axes):
-        network_builder.cells[gid].plot_voltage(ax)
-        ax.set_title('%s (gid=%d)' % (cells[gid].name, gid))
+# Let's do this for the rest of the cell types
+fig, axes = plt.subplots(1, 2, sharey=True, figsize=(8, 4))
+for gid, ax in zip([gid_dict['L2_pyramidal'][0],
+                    gid_dict['L5_pyramidal'][0]], axes):
+    ax.plot(times, vsoma[gid])
+    ax.set_title('%s (gid=%d)' % (net.gid_to_type(gid), gid))
+    ax.set_xlabel('Time (ms)')
+    ax.set_ylabel('Voltage (mV)')
+plt.show()
+
+###############################################################################
+# The spiking activity across cell types can also be visualized with raster
+# plots and histograms
+
+# Spike raster plot
+net.spikes.plot()
+
+# Spike histogram
+fig, axes = plt.subplots(1, 1, figsize=(8, 4))
+net.spikes.plot_hist(ax=axes, spike_types=['L2_pyramidal', 'L5_pyramidal'])
+axes.set_xlabel('Time (ms)')

hnn_core/cell.py

@@ -67,6 +67,9 @@ class _Cell(ABC):
         The Dipole objects (see dipole.mod).
     dict_currents : dict of h.Vector()
         The soma currents (keys are soma_gabaa, soma_gabab etc.)
+    rec_v : h.Vector()
+        Recording of somatic voltage. Must be enabled
+        by running simulate_dipole(net, record_vsoma=True)
     """
 
     def __init__(self, gid, soma_props):
@@ -75,6 +78,7 @@ def __init__(self, gid, soma_props):
         self.list_IClamp = None
         self.soma_props = soma_props
         self.create_soma()
+        self.rec_v = h.Vector()
 
     def __repr__(self):
         class_name = self.__class__.__name__
@@ -207,6 +211,10 @@ def record_current_soma(self):
                 " but no self.synapses dict was found")
             pass
 
+    def record_voltage_soma(self):
+        """Record voltage at soma."""
+        self.rec_v.record(self.soma(0.5)._ref_v)
+
     def syn_create(self, secloc, e, tau1, tau2):
         """Create an h.Exp2Syn synapse.
 
@@ -298,62 +306,3 @@ def shape_soma(self):
         # self.soma.diam set above
         h.pt3dadd(0, 0, 0, self.diam, sec=self.soma)
         h.pt3dadd(0, self.L, 0, self.diam, sec=self.soma)
-
-    def plot_voltage(self, ax=None, delay=2, duration=100, dt=0.2,
-                     amplitude=1, show=True):
-        """Plot voltage on soma for an injected current
-
-        Parameters
-        ----------
-        ax : instance of matplotlib axis | None
-            An axis object from matplotlib. If None,
-            a new figure is created.
-        delay : float (in ms)
-            The start time of the injection current.
-        duration : float (in ms)
-            The duration of the injection current.
-        dt : float (in ms)
-            The integration time step
-        amplitude : float (in nA)
-            The amplitude of the injection current.
-        show : bool
-            Call plt.show() if True. Set to False if working in
-            headless mode (e.g., over a remote connection).
-        """
-        import matplotlib.pyplot as plt
-        from neuron import h
-        h.load_file('stdrun.hoc')
-
-        soma = self.soma
-
-        h.tstop = duration
-        h.dt = dt
-        h.celsius = 37
-
-        iclamp = h.IClamp(soma(0.5))
-        iclamp.delay = delay
-        iclamp.dur = duration
-        iclamp.amp = amplitude
-
-        v_membrane = h.Vector().record(self.soma(0.5)._ref_v)
-        times = h.Vector().record(h._ref_t)
-
-        print('Simulating soma voltage')
-        h.finitialize()
-
-        def simulation_time():
-            print('Simulation time: {0} ms...'.format(round(h.t, 2)))
-
-        for tt in range(0, int(h.tstop), 10):
-            h.CVode().event(tt, simulation_time)
-
-        h.continuerun(duration)
-        print('[Done]')
-
-        if ax is None:
-            fig, ax = plt.subplots(1, 1)
-        ax.plot(times, v_membrane)
-        ax.set_xlabel('Time (ms)')
-        ax.set_ylabel('Voltage (mV)')
-        if show:
-            plt.show()

hnn_core/dipole.py

@@ -17,7 +17,7 @@ def _hammfilt(x, winsz):
     return convolve(x, win, 'same')
 
 
-def simulate_dipole(net, n_trials=None):
+def simulate_dipole(net, n_trials=None, record_vsoma=False):
     """Simulate a dipole given the experiment parameters.
 
     Parameters
@@ -28,6 +28,8 @@ def simulate_dipole(net, n_trials=None):
     n_trials : int | None
         The number of trials to simulate. If None the value in
         net.params['N_trials'] will be used
+    record_vsoma : bool
+        Option to record somatic voltages from cells
 
     Returns
     -------
@@ -48,6 +50,12 @@ def simulate_dipole(net, n_trials=None):
     if n_trials < 1:
         raise ValueError("Invalid number of simulations: %d" % n_trials)
 
+    if record_vsoma is not None and isinstance(record_vsoma, bool):
+        net.params['record_vsoma'] = record_vsoma
+    else:
+        raise TypeError("record_vsoma must be bool, got %s"
+                        % type(record_vsoma).__name__)
+
     dpls = _BACKEND.simulate(net)
 
     return dpls