Record calcium concentration new

doc/whats_new.rst

@@ -200,7 +200,7 @@ API
   :meth:`~hnn_core.extracellular.ExtracellularArray.plot_lfp`, by
   `Steven Brandt`_ and `Ryan Thorpe`_ in :gh:`517`.
 
-- Recorded voltages/currents from the soma, as well all sections, are enabled by
+- Recorded voltages/currents from the soma, as well as all sections, are enabled by
   setting either `record_vsec` or `record_isec` to `'all'` or `'soma'` 
   in :func:`~hnn_core.simulate_dipole`. Recordings are now accessed through
   :class:`~hnn_core.CellResponse.vsec` and :class:`~hnn_core.CellResponse.isec`,
@@ -210,6 +210,11 @@ API
   :class:`~hnn_core.Network` objects, 
   by `Nick Tolley`_ and `Ryan Thorpe`_ in :gh:`619`.
 
+- Recorded calcium conncetration from the soma, as well as all sections, are enabled
+  by setting `record_ca` to `soma` or `all` in :func:`~hnn_core.simulate_dipole`.
+  Recordings are accessed through :class:`~hnn_core.CellResponse.ca`, 
+  by `Katharina Duecker`_ in :gh:`792`
+
 People who contributed to this release (in alphabetical order):
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

hnn_core/cell.py

@@ -348,6 +348,9 @@ class Cell:
         by synapse type (keys can be soma_gabaa, soma_gabab etc.).
         Must be enabled by running simulate_dipole(net, record_isec=True)
         or simulate_dipole(net, record_isoma=True)
+    ca : dict
+        Contains recording of section speicifc calcium concentration.
+        Must be enabled by running simulate_dipole(net, record_ca=True).
     tonic_biases : list of h.IClamp
         The current clamps inserted at each section of the cell
         for tonic biasing inputs.
@@ -391,6 +394,7 @@ def __init__(self, name, pos, sections, synapses, sect_loc, cell_tree,
         self.dipole_pp = list()
         self.vsec = dict()
         self.isec = dict()
+        self.ca = dict()
         # insert iclamp
         self.list_IClamp = list()
         self._gid = None
@@ -468,6 +472,7 @@ def to_dict(self):
         cell_data['dipole_pp'] = self.dipole_pp
         cell_data['vsec'] = self.vsec
         cell_data['isec'] = self.isec
+        cell_data['ca'] = self.ca
         cell_data['tonic_biases'] = self.tonic_biases
         return cell_data
 
@@ -743,7 +748,7 @@ def create_tonic_bias(self, amplitude, t0, tstop, loc=0.5):
         stim.amp = amplitude
         self.tonic_biases.append(stim)
 
-    def record(self, record_vsec=False, record_isec=False):
+    def record(self, record_vsec=False, record_isec=False, record_ca=False):
         """ Record current and voltage from all sections
 
         Parameters
@@ -754,6 +759,9 @@ def record(self, record_vsec=False, record_isec=False):
         record_isec : 'all' | 'soma' | False
             Option to record voltages from all sections ('all'), or just
             the soma ('soma'). Default: False.
+        record_ca : 'all' | 'soma' | False
+            Option to record calcium concentration from all sections ('all'),
+            or just the soma ('soma'). Default: False.
         """
 
         section_names = list(self.sections.keys())
@@ -783,10 +791,35 @@ def record(self, record_vsec=False, record_isec=False):
 
                 for syn_name in self.isec[sec_name]:
                     self.isec[sec_name][syn_name] = h.Vector()
-
                     self.isec[sec_name][syn_name].record(
                         self._nrn_synapses[syn_name]._ref_i)
 
+        # calcium concentration
+        if record_ca == 'soma':
+            self.ca = dict.fromkeys(['soma'])
+        elif record_ca == 'all':
+            self.ca = dict.fromkeys(section_names)
+
+        if record_ca:
+            for sec_name in self.ca:
+                if hasattr(self._nrn_sections[sec_name](0.5), '_ref_cai'):
+                    self.ca[sec_name] = h.Vector()
+                    self.ca[sec_name].record(
+                        self._nrn_sections[sec_name](0.5)._ref_cai)
+
+        # calcium concentration
+        if record_ca == 'soma':
+            self.ca = dict.fromkeys(['soma'])
+        elif record_ca == 'all':
+            self.ca = dict.fromkeys(section_names)
+
+        if record_ca:
+            for sec_name in self.ca:
+                if hasattr(self._nrn_sections[sec_name](0.5), '_ref_cai'):
+                    self.ca[sec_name] = h.Vector()
+                    self.ca[sec_name].record(
+                        self._nrn_sections[sec_name](0.5)._ref_cai)
+
     def syn_create(self, secloc, e, tau1, tau2):
         """Create an h.Exp2Syn synapse.
 

hnn_core/cell_response.py

@@ -56,6 +56,10 @@ class CellResponse(object):
     isec : list (n_trials,) of dict, shape
         Each element of the outer list is a trial.
         Dictionary indexed by gids containing currents for cell sections.
+    ca : list (n_trials,) of dict, shape
+        Each element of the outer list is a trial.
+        Dictionary indexed by gids containing calcium concentration
+        for cell sections.
     times : array-like, shape (n_times,)
         Array of time points for samples in continuous data.
         This includes vsoma and isoma.
@@ -115,6 +119,7 @@ def __init__(self, spike_times=None, spike_gids=None, spike_types=None,
         self._spike_types = spike_types
         self._vsec = list()
         self._isec = list()
+        self._ca = list()
         if times is not None:
             if not isinstance(times, (list, np.ndarray)):
                 raise TypeError("'times' is an np.ndarray of simulation times")
@@ -138,7 +143,8 @@ def __eq__(self, other):
                 self._spike_gids == other._spike_gids and
                 self._spike_types == other._spike_types and
                 self.vsec == other.vsec and
-                self.isec == other.isec)
+                self.isec == other.isec and
+                self.ca == other.ca)
 
     def __getitem__(self, gid_item):
         """Returns a CellResponse object with a copied subset filtered by gid.
@@ -228,6 +234,10 @@ def vsec(self):
     def isec(self):
         return self._isec
 
+    @property
+    def ca(self):
+        return self._ca
+
     @property
     def times(self):
         return self._times
@@ -427,6 +437,12 @@ def to_dict(self):
             # Turn `int` gid keys into string values for hdf5 format
             trial = dict((str(key), val) for key, val in trial.items())
             cell_response_data['isec'].append(trial)
+        ca_data = self.ca
+        cell_response_data['ca'] = list()
+        for trial in ca_data:
+            # Turn `int` gid keys into string values for hdf5 format
+            trial = dict((str(key), val) for key, val in trial.items())
+            cell_response_data['ca'].append(trial)
         cell_response_data['times'] = self.times
         return cell_response_data
 

hnn_core/dipole.py

@@ -16,7 +16,7 @@
 
 
 def simulate_dipole(net, tstop, dt=0.025, n_trials=None, record_vsec=False,
-                    record_isec=False, postproc=False):
+                    record_isec=False, record_ca=False, postproc=False):
     """Simulate a dipole given the experiment parameters.
 
     Parameters
@@ -35,8 +35,11 @@ def simulate_dipole(net, tstop, dt=0.025, n_trials=None, record_vsec=False,
         Option to record voltages from all sections ('all'), or just
         the soma ('soma'). Default: False.
     record_isec : 'all' | 'soma' | False
-        Option to record voltages from all sections ('all'), or just
+        Option to record synaptic currents from all sections ('all'), or just
         the soma ('soma'). Default: False.
+    record_ca : 'all' | 'soma' | False
+        Option to record calcium concentration from all sections ('all'),
+        or just the soma ('soma'). Default: False.
     postproc : bool
         If True, smoothing (``dipole_smooth_win``) and scaling
         (``dipole_scalefctr``) values are read from the parameter file, and
@@ -96,6 +99,10 @@ def simulate_dipole(net, tstop, dt=0.025, n_trials=None, record_vsec=False,
 
     net._params['record_isec'] = record_isec
 
+    _check_option('record_ca', record_ca, ['all', 'soma', False])
+
+    net._params['record_ca'] = record_ca
+
     if postproc:
         warnings.warn('The postproc-argument is deprecated and will be removed'
                       ' in a future release of hnn-core. Please define '

hnn_core/network_builder.py

@@ -100,6 +100,13 @@ def simulation_time():
             isec_py[gid][sec_name] = {
                 key: isec.to_python() for key, isec in isec.items()}
 
+    ca_py = dict()
+    for gid, ca_dict in neuron_net._ca.items():
+        ca_py[gid] = dict()
+        for sec_name, ca in ca_dict.items():
+            if ca is not None:
+                ca_py[gid][sec_name] = ca.to_python()
+
     dpl_data = np.c_[
         neuron_net._nrn_dipoles['L2_pyramidal'].as_numpy() +
         neuron_net._nrn_dipoles['L5_pyramidal'].as_numpy(),
@@ -119,6 +126,7 @@ def simulation_time():
             'gid_ranges': net.gid_ranges,
             'vsec': vsec_py,
             'isec': isec_py,
+            'ca': ca_py,
             'rec_data': rec_arr_py,
             'rec_times': rec_times_py,
             'times': times.to_python()}
@@ -291,6 +299,7 @@ def __init__(self, net, trial_idx=0):
 
         self._vsec = dict()
         self._isec = dict()
+        self._ca = dict()
         self._nrn_rec_arrays = dict()
         self._nrn_rec_callbacks = list()
 
@@ -327,9 +336,11 @@ def _build(self):
 
         record_vsec = self.net._params['record_vsec']
         record_isec = self.net._params['record_isec']
+        record_ca = self.net._params['record_ca']
         self._create_cells_and_drives(threshold=self.net._params['threshold'],
                                       record_vsec=record_vsec,
-                                      record_isec=record_isec)
+                                      record_isec=record_isec,
+                                      record_ca=record_ca)
 
         self.state_init()
 
@@ -399,7 +410,7 @@ def _gid_assign(self, rank=None, n_hosts=None):
         self._gid_list.sort()
 
     def _create_cells_and_drives(self, threshold, record_vsec=False,
-                                 record_isec=False):
+                                 record_isec=False, record_ca=False):
         """Parallel create cells AND external drives
 
         NB: _Cell.__init__ calls h.Section -> non-picklable!
@@ -434,7 +445,7 @@ def _create_cells_and_drives(self, threshold, record_vsec=False,
                         src_type in self.net.external_biases['tonic']):
                     cell.create_tonic_bias(**self.net.external_biases
                                            ['tonic'][src_type])
-                cell.record(record_vsec, record_isec)
+                cell.record(record_vsec, record_isec, record_ca)
 
                 # this call could belong in init of a _Cell (with threshold)?
                 nrn_netcon = cell.setup_source_netcon(threshold)
@@ -564,6 +575,7 @@ def aggregate_data(self, n_samples):
 
             self._vsec[cell.gid] = cell.vsec
             self._isec[cell.gid] = cell.isec
+            self._ca[cell.gid] = cell.ca
 
         # reduce across threads
         for nrn_dpl in self._nrn_dipoles.values():
@@ -574,6 +586,7 @@ def aggregate_data(self, n_samples):
         # aggregate the currents and voltages independently on each proc
         vsec_list = _PC.py_gather(self._vsec, 0)
         isec_list = _PC.py_gather(self._isec, 0)
+        ca_list = _PC.py_gather(self._ca, 0)
 
         # combine spiking data from each proc
         spike_times_list = _PC.py_gather(self._spike_times, 0)
@@ -589,6 +602,8 @@ def aggregate_data(self, n_samples):
                 self._vsec.update(vsec)
             for isec in isec_list:
                 self._isec.update(isec)
+            for ca in ca_list:
+                self._ca.update(ca)
 
         _PC.barrier()  # get all nodes to this place before continuing
 

hnn_core/parallel_backends.py

@@ -53,6 +53,7 @@ def _gather_trial_data(sim_data, net, n_trials, postproc):
         net.cell_response.update_types(net.gid_ranges)
         net.cell_response._vsec.append(sim_data[idx]['vsec'])
         net.cell_response._isec.append(sim_data[idx]['isec'])
+        net.cell_response._ca.append(sim_data[idx]['ca'])
 
         # extracellular array
         for arr_name, arr in net.rec_arrays.items():

hnn_core/params_default.py

@@ -149,6 +149,7 @@ def get_params_default(nprox=2, ndist=1):
         'record_isoma': 0,  # whether to record somatic currents
         'record_vsec': 0,  # whether to record voltages
         'record_isec': 0,  # whether to record currents
+        'record_ca': 0,  # whether to record calcium concentration
 
         # numerics
         # N_trials of 1 means that seed is set by rank

hnn_core/tests/conftest.py

@@ -75,8 +75,8 @@ def pytest_runtest_setup(item):
 def run_hnn_core_fixture():
     def _run_hnn_core_fixture(backend=None, n_procs=None, n_jobs=1,
                               reduced=False, record_vsec=False,
-                              record_isec=False, postproc=False,
-                              electrode_array=None):
+                              record_isec=False, record_ca=False,
+                              postproc=False, electrode_array=None):
         hnn_core_root = op.dirname(hnn_core.__file__)
 
         # default params
@@ -106,15 +106,18 @@ def _run_hnn_core_fixture(backend=None, n_procs=None, n_jobs=1,
             with MPIBackend(n_procs=n_procs, mpi_cmd='mpiexec'):
                 dpls = simulate_dipole(net, record_vsec=record_vsec,
                                        record_isec=record_isec,
+                                       record_ca=record_ca,
                                        postproc=postproc, tstop=tstop)
         elif backend == 'joblib':
             with JoblibBackend(n_jobs=n_jobs):
                 dpls = simulate_dipole(net, record_vsec=record_vsec,
                                        record_isec=record_isec,
+                                       record_ca=record_ca,
                                        postproc=postproc, tstop=tstop)
         else:
             dpls = simulate_dipole(net, record_vsec=record_vsec,
                                    record_isec=record_isec,
+                                   record_ca=record_ca,
                                    postproc=postproc, tstop=tstop)
 
         # check that the network object is picklable after the simulation

hnn_core/tests/test_cell_response.py

@@ -47,7 +47,7 @@ def test_cell_response(tmp_path):
     # reset clears all recorded variables, but leaves simulation time intact
     assert len(cell_response.times) == len(sim_times)
     sim_attributes = ['_spike_times', '_spike_gids', '_spike_types',
-                      '_vsec', '_isec']
+                      '_vsec', '_isec', '_ca']
     net_attributes = ['_times', '_cell_type_names']  # `Network.__init__`
     # creates these check that we always know which response attributes are
     # simulated see #291 for discussion; objective is to keep cell_response

hnn_core/tests/test_dipole.py

@@ -151,13 +151,14 @@ def test_dipole(tmp_path, run_hnn_core_fixture):
     # XXX all below to be deprecated in 0.3
     dpls_raw, net = run_hnn_core_fixture(backend='joblib', n_jobs=1,
                                          reduced=True, record_isec='soma',
-                                         record_vsec='soma')
+                                         record_vsec='soma', record_ca='soma')
     # test deprecation of postproc
     with pytest.warns(DeprecationWarning,
                       match='The postproc-argument is deprecated'):
         dpls, _ = run_hnn_core_fixture(backend='joblib', n_jobs=1,
                                        reduced=True, record_isec='soma',
-                                       record_vsec='soma', postproc=True)
+                                       record_vsec='soma', record_ca='soma',
+                                       postproc=True)
     with pytest.raises(AssertionError):
         assert_allclose(dpls[0].data['agg'], dpls_raw[0].data['agg'])
 
@@ -186,6 +187,9 @@ def test_dipole_simulation():
     with pytest.raises(ValueError, match="Invalid value for the"):
         simulate_dipole(net, tstop=25., n_trials=1, record_vsec=False,
                         record_isec='abc')
+    with pytest.raises(ValueError, match="Invalid value for the"):
+        simulate_dipole(net, tstop=25., n_trials=1, record_vsec=False,
+                        record_isec=False, record_ca='abc')
 
     # test Network.copy() returns 'bare' network after simulating
     dpl = simulate_dipole(net, tstop=25., n_trials=1)[0]
@@ -214,20 +218,22 @@ def test_cell_response_backends(run_hnn_core_fixture):
     trial_idx, n_trials, gid = 0, 2, 7
     _, joblib_net = run_hnn_core_fixture(backend='joblib', n_jobs=1,
                                          reduced=True, record_vsec='all',
-                                         record_isec='soma')
+                                         record_isec='soma', record_ca='all')
     _, mpi_net = run_hnn_core_fixture(backend='mpi', n_procs=2, reduced=True,
-                                      record_vsec='all', record_isec='soma')
+                                      record_vsec='all', record_isec='soma',
+                                      record_ca='all')
+
     n_times = len(joblib_net.cell_response.times)
 
     assert len(joblib_net.cell_response.vsec) == n_trials
     assert len(joblib_net.cell_response.isec) == n_trials
+    assert len(joblib_net.cell_response.ca) == n_trials
     assert len(joblib_net.cell_response.vsec[trial_idx][gid]) == 8  # num sec
     assert len(joblib_net.cell_response.isec[trial_idx][gid]) == 1
     assert len(joblib_net.cell_response.vsec[
         trial_idx][gid]['apical_1']) == n_times
     assert len(joblib_net.cell_response.isec[
                trial_idx][gid]['soma']['soma_gabaa']) == n_times
-
     assert len(mpi_net.cell_response.vsec) == n_trials
     assert len(mpi_net.cell_response.isec) == n_trials
     assert len(mpi_net.cell_response.vsec[trial_idx][gid]) == 8  # num sec
@@ -239,6 +245,16 @@ def test_cell_response_backends(run_hnn_core_fixture):
     assert mpi_net.cell_response.vsec == joblib_net.cell_response.vsec
     assert mpi_net.cell_response.isec == joblib_net.cell_response.isec
 
+    # test if calcium concentration is stored correctly (only L5 pyramidal)
+    gid = joblib_net.gid_ranges['L5_pyramidal'][0]
+    assert len(joblib_net.cell_response.ca[trial_idx][gid]) == 9
+    assert len(joblib_net.cell_response.ca[
+               trial_idx][gid]['soma']) == n_times
+    assert len(mpi_net.cell_response.ca[trial_idx][gid]) == 9
+    assert len(mpi_net.cell_response.ca[
+               trial_idx][gid]['soma']) == n_times
+    assert len(mpi_net.cell_response.ca[
+               trial_idx][gid]['apical_1']) == n_times
     # Test if spike time falls within depolarization window above v_thresh
     v_thresh = 0.0
     times = np.array(joblib_net.cell_response.times)