[MRG] ENH: calcium dynamics

doc/api.rst

@@ -27,10 +27,18 @@ Published Models (:py:mod:`hnn_core`):
 .. autosummary::
    :toctree: generated/
 
-   default_network
    jones_2009_model
    law_2021_model
 
+Unpublished Models (:py:mod:`hnn_core`):
+----------------------------------------
+.. currentmodule:: hnn_core
+
+.. autosummary::
+   :toctree: generated/
+
+   calcium_model
+
 Dipole (:py:mod:`hnn_core.dipole`):
 -----------------------------------
 

doc/whats_new.rst

@@ -34,6 +34,8 @@ Changelog
 
 - Add ability to interactivity explore connections in :func:`~hnn_core.viz.plot_cell_connectivity` by `Mainak Jas`_ in `#376 <https://github.com/jonescompneurolab/hnn-core/pull/376>`_
 
+- Add ``calcium_model`` with a distance dependent calcium channel conductivity, by `Nick Tolley`_ in `#348 <https://github.com/jonescompneurolab/hnn-core/pull/333>`_
+
 Bug
 ~~~
 
@@ -59,7 +61,7 @@ API
 - 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>`_
 
-- New API for network creation. The default network is now created with ``net = default_network(params)``, by `Nick Tolley`_ in `#318 <https://github.com/jonescompneurolab/hnn-core/pull/318>`_
+- New API for network creation. The default network is now created with ``net = jones_2009_model(params)``, by `Nick Tolley`_ in `#318 <https://github.com/jonescompneurolab/hnn-core/pull/318>`_
 
 - Replace parameter `T` with `tstop` in :func:`~hnn_core.Network.add_tonic_bias` and :func:`~hnn_core.Cell.create_tonic_bias` to be more consistent with other functions and improve readability, by `Kenneth Loi`_ in `#354 <https://github.com/jonescompneurolab/hnn-core/pull/354>`_
 

examples/howto/plot_connectivity.py

@@ -16,7 +16,7 @@
 # Let us import ``hnn_core``.
 
 import hnn_core
-from hnn_core import read_params, default_network, simulate_dipole
+from hnn_core import read_params, jones_2009_model, simulate_dipole
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 
@@ -32,7 +32,7 @@
 # explore how it changes with new connections. We first instantiate the
 # network. (Note: Setting ``add_drives_from_params=True`` loads a set of
 # predefined drives without the drives API shown previously).
-net_erp = default_network(params, add_drives_from_params=True)
+net_erp = jones_2009_model(params, add_drives_from_params=True)
 
 ###############################################################################
 # Instantiating the network comes with a predefined set of connections that
@@ -45,7 +45,7 @@
 
 print(len(net_erp.connectivity))
 
-conn_idx = 20
+conn_idx = 6
 print(net_erp.connectivity[conn_idx])
 plot_connectivity_matrix(net_erp, conn_idx)
 
@@ -73,7 +73,7 @@
 
 
 def get_network(probability=1.0):
-    net = default_network(params, add_drives_from_params=True)
+    net = jones_2009_model(params, add_drives_from_params=True)
     net.clear_connectivity()
 
     # Pyramidal cell connections

examples/howto/plot_firing_pattern.py

@@ -17,7 +17,8 @@
 # Let us import ``hnn_core``.
 
 import hnn_core
-from hnn_core import read_params, read_spikes, default_network, simulate_dipole
+from hnn_core import (read_params, read_spikes, jones_2009_model,
+                      simulate_dipole)
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 
@@ -31,7 +32,7 @@
 # :ref:`evoked example <sphx_glr_auto_examples_plot_simulate_evoked.py>`.
 import matplotlib.pyplot as plt
 
-net = default_network(params)
+net = jones_2009_model(params)
 
 ###############################################################################
 # ``net`` does not have any driving inputs and only defines the local network

examples/howto/plot_record_extracellular_potentials.py

@@ -29,12 +29,12 @@
 # file.
 
 import hnn_core
-from hnn_core import read_params, default_network, simulate_dipole
+from hnn_core import read_params, jones_2009_model, simulate_dipole
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 params_fname = op.join(hnn_core_root, 'param', 'default.json')
 params = read_params(params_fname)
-net = default_network(params, add_drives_from_params=True)
+net = jones_2009_model(params, add_drives_from_params=True)
 
 ###############################################################################
 # Extracellular recordings require specifying the electrode postions. It can be

examples/howto/plot_simulate_mpi_backend.py

@@ -21,7 +21,7 @@
 import os.path as op
 
 import hnn_core
-from hnn_core import simulate_dipole, read_params, default_network
+from hnn_core import simulate_dipole, read_params, jones_2009_model
 
 ###############################################################################
 # Then we setup the directories and Neuron
@@ -47,7 +47,7 @@
 # The occurrence of each burst is jittered by a random, normally distributed
 # amount (20 ms standard deviation). We repeat the burst train 10 times, each
 # time with unique randomization.
-net = default_network(params)
+net = jones_2009_model(params)
 
 weights_ampa = {'L2_pyramidal': 5.4e-5, 'L5_pyramidal': 5.4e-5}
 net.add_bursty_drive(

examples/workflows/plot_simulate_alpha.py

@@ -26,7 +26,7 @@
 # Let us import hnn_core
 
 import hnn_core
-from hnn_core import simulate_dipole, read_params, default_network
+from hnn_core import simulate_dipole, read_params, jones_2009_model
 
 ###############################################################################
 # Then we setup the directories and read the default parameters file
@@ -43,7 +43,7 @@
 # time with unique randomization. The drive is only connected to the proximal
 # (dendritic) AMPA synapses on L2/3 and L5 pyramidal neurons.
 params['tstop'] = 310
-net = default_network(params)
+net = jones_2009_model(params)
 
 location = 'proximal'
 burst_std = 20

examples/workflows/plot_simulate_beta.py

@@ -1,14 +1,15 @@
 """
 ===============================
-06. Simulate beta modulated ERP
+05. Simulate beta modulated ERP
 ===============================
 
 This example demonstrates how event related potentials (ERP) are modulated
-by prestimulus beta events. Transient beta activity in the neocortex has
-been shown to suppress the perceptibility of sensory input. This suppression
-depends on the timing of the beta event, and the incoming sensory
-information. The following example demonstrates the biophysical mechansisms
-underlying beta mediated sensory suppression.
+by prestimulus beta events. Specifically, this example reproduces Figure 5
+from Law et al. 2021 [1]_. To be consistent with the publication, the default
+network connectivity is altered. These modfications demonstrate a potential
+mechanism by which transient beta activity in the neocortex can suppress
+the perceptibility of sensory input. This suppression depends on the timing
+of the beta event, and the incoming sensory information.
 """
 
 # Authors: Nick Tolley <[email protected]>
@@ -48,7 +49,7 @@
 # A major change to the Jones 2009 model is the addition of a
 # Martinotti-like recurrent tuft connection [3]_. This new connection
 # originates from L5 basket cells, and provides GABAa inhibition on
-# the distal dendrites of L5 basket cells.
+# the distal dendrites of L5 pyramidal cells.
 print('Recurrent Tuft Connection')
 print(net.connectivity[16])
 
@@ -69,12 +70,13 @@
 # but modified to reflect the parameters used in Law et al. 2021.
 # Specifically, we are considering the case where a tactile stimulus is
 # delivered at 150 ms. 25 ms later, the first input to sensory cortex arrives
-# as at the proximal drive to the cortical column. Proximal drive corresponds
-# projects from the direct thalamic nuclei. This is followed by one distal
-# representing projections from indirect thalamic nuclei, and a final late
-# proximal drive. It is important to note that the parameter values for each
-# are different from previous examples of the evoekd response. This reflects
-# the altered network dynamics due to the changes described above.
+# as a proximal drive to the cortical column. Proximal drive corresponds to
+# projections from the direct thalamic nuclei. This is followed by one distal
+# drive representing projections from indirect thalamic nuclei, and a final
+# late proximal drive. It is important to note that the parameter values for
+# each are different from previous examples of the evoked response.
+# This reflects the altered network dynamics due to the changes described
+# above.
 def add_erp_drives(net, stimulus_start):
     # Distal evoked drive
     weights_ampa_d1 = {'L2_basket': 0.0005, 'L2_pyramidal': 0.004,
@@ -181,15 +183,15 @@ def add_beta_drives(net, beta_start):
 net_beta.cell_response.plot_spikes_raster(ax=axes[2], show=False)
 axes[2].set_title('Spike Raster')
 
-# Create an fixed-step tiling of frequencies from 1 to 40 Hz in steps of 1 Hz
+# Create a fixed-step tiling of frequencies from 1 to 40 Hz in steps of 1 Hz
 freqs = np.arange(10., 60., 1.)
 dpls_beta[0].plot_tfr_morlet(freqs, n_cycles=7, ax=axes[3])
 
 ###############################################################################
 # Next we will inspect what happens when a sensory stimulus is delivered 75 ms
 # after a beta event. Note that the delay time for a tactile stimulus at the
 # hand to arrive at the cortex is roughly 25 ms, which means the first proximal
-# input to thecortical column occurs ~100 ms after the beta event.
+# input to the cortical column occurs ~100 ms after the beta event.
 dpls_beta_erp[0].smooth(45)
 fig, axes = plt.subplots(3, 1, sharex=True, figsize=(7, 7),
                          constrained_layout=True)
@@ -201,11 +203,11 @@ def add_beta_drives(net, beta_start):
 axes[2].set_title('Spike Raster')
 
 ###############################################################################
-# To help understand the effect of beta mediated inhibition of the response to
+# To help understand the effect of beta mediated inhibition on the response to
 # incoming sensory stimuli, we can compare the ERP and spiking activity due to
 # sensory input with and without a beta event.
-# The sustained inhibition of the network ultimately depressing
-# the sensory response which is assoicated with a reduced ERP amplitude
+# The sustained inhibition of the network ultimately depresses
+# the sensory response which is associated with a reduced ERP amplitude
 dpls_erp[0].smooth(45)
 fig, axes = plt.subplots(3, 1, sharex=True, figsize=(7, 7),
                          constrained_layout=True)

examples/workflows/plot_simulate_evoked.py

@@ -30,7 +30,7 @@
 # Let us import hnn_core
 
 import hnn_core
-from hnn_core import simulate_dipole, read_params, default_network
+from hnn_core import simulate_dipole, read_params, jones_2009_model
 from hnn_core.viz import plot_dipole
 
 hnn_core_root = op.dirname(hnn_core.__file__)
@@ -50,7 +50,7 @@
 ###############################################################################
 # Let us first create our network from the params file and visualize the cells
 # inside it.
-net = default_network(params)
+net = jones_2009_model(params)
 net.plot_cells()
 net.cell_types['L5_pyramidal'].plot_morphology()
 

examples/workflows/plot_simulate_gamma.py

@@ -27,7 +27,7 @@
 # Let us import hnn_core
 
 import hnn_core
-from hnn_core import simulate_dipole, read_params, default_network
+from hnn_core import simulate_dipole, read_params, jones_2009_model
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 
@@ -44,7 +44,7 @@
 # simulate the dipole moment in a single trial (the default value used by
 # ``simulate_dipole`` is ``n_trials=params['N_trials']``).
 
-net = default_network(params)
+net = jones_2009_model(params)
 
 weights_ampa = {'L2_pyramidal': 0.0008, 'L5_pyramidal': 0.0075}
 synaptic_delays = {'L2_pyramidal': 0.1, 'L5_pyramidal': 1.0}

examples/workflows/plot_simulate_somato.py

@@ -143,15 +143,15 @@
 # network parameters from ``N20.json`` and instantiate the network.
 
 import hnn_core
-from hnn_core import simulate_dipole, read_params, default_network
+from hnn_core import simulate_dipole, read_params, jones_2009_model
 from hnn_core import average_dipoles, JoblibBackend
 
 hnn_core_root = op.dirname(hnn_core.__file__)
 
 params_fname = op.join(hnn_core_root, 'param', 'N20.json')
 params = read_params(params_fname)
 
-net = default_network(params)
+net = jones_2009_model(params)
 
 ###############################################################################
 # To simulate the source of the median nerve evoked response, we add a

hnn_core/__init__.py

@@ -2,7 +2,7 @@
 from .drives import drive_event_times
 from .params import Params, read_params
 from .network import Network
-from .network_models import default_network, jones_2009_model, law_2021_model
+from .network_models import jones_2009_model, law_2021_model, calcium_model
 from .cell import Cell
 from .cell_response import CellResponse, read_spikes
 from .cells_default import pyramidal, basket

hnn_core/cells_default.py

@@ -4,7 +4,7 @@
 #          Sam Neymotin <[email protected]>
 
 import numpy as np
-
+from functools import partial
 from .cell import Cell
 
 from .params import compare_dictionaries
@@ -280,3 +280,36 @@ def pyramidal(cell_name, pos=(0, 0, 0), override_params=None, gid=None):
                 topology=topology,
                 sect_loc=sect_loc,
                 gid=gid)
+
+
+def _get_g_at_dist(x, gsoma, gdend, xkink):
+    """Compute distance-dependent ionic conductance.
+
+    Notes
+    -----
+    Linearly scales conductance along dendrite.
+    Returns gdend when x > xkink.
+    """
+    return gsoma + np.min([xkink, x]) * (gdend - gsoma) / xkink
+
+
+def pyramidal_ca(cell_name, pos, override_params=None, gid=None):
+    """Calcium dynamics."""
+
+    if override_params is None:
+        override_params = dict()
+
+    override_params['L5Pyr_soma_gkbar_hh2'] = 0.06
+    override_params['L5Pyr_soma_gnabar_hh2'] = 0.32
+    override_params['L5Pyr_dend_gkbar_hh2'] = 1e-4
+    override_params['L5Pyr_dend_gnabar_hh2'] = 28e-4
+
+    gbar_ca = partial(_get_g_at_dist,
+                      gsoma=10.,
+                      gdend=40.,
+                      xkink=1501)
+    override_params['L5Pyr_dend_gbar_ca'] = gbar_ca
+    cell = pyramidal(cell_name, pos, override_params=override_params,
+                     gid=gid)
+
+    return cell

hnn_core/network.py

@@ -211,7 +211,7 @@ class Network(object):
 
     Notes
     ----
-    `net = default_network(params)` is the reccomended path for creating a
+    `net = jones_2009_model(params)` is the reccomended path for creating a
     network. Instantiating the network as `net = Network(params)` will
     produce a network with no cell to cell connections. As such,
     connectivity information contained in `params` will be ignored.