.. currentmodule:: brian2

.. Destexhe_et_al_1998:

Example: Destexhe_et_al_1998
============================


        .. only:: html

            .. |launchbinder| image:: file:///usr/share/doc/python-brian-doc/docs/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/frompapers/Destexhe_et_al_1998.ipynb

            .. note::
               You can launch an interactive, editable version of this
               example without installing any local files
               using the Binder service (although note that at some times this
               may be slow or fail to open): |launchbinder|_

        

Reproduces Figure 12 (simplified three-compartment model) from the following
paper:

Dendritic Low-Threshold Calcium Currents in Thalamic Relay Cells
Alain Destexhe, Mike Neubig, Daniel Ulrich, John Huguenard
Journal of Neuroscience 15 May 1998, 18 (10) 3574-3588

The original NEURON code is available on ModelDB: https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=279

Reference for the original morphology:

Rat VB neuron (thalamocortical cell), given by J. Huguenard, stained with
biocytin and traced by A. Destexhe, December 1992.  The neuron is described
in: J.R. Huguenard & D.A. Prince, A novel T-type current underlies prolonged
calcium-dependent burst firing in GABAergic neurons of rat thalamic reticular
nucleus.  J. Neurosci. 12: 3804-3817, 1992.

Available at NeuroMorpho.org:

http://neuromorpho.org/neuron_info.jsp?neuron_name=tc200
NeuroMorpho.Org ID :NMO_00881

Notes
-----
* Completely removed the "Fast mechanism for submembranal Ca++ concentration
  (cai)" -- it did not affect the results presented here
* Time constants for the I_T current are slightly different from the equations
  given in the paper -- the paper calculation seems to be based on 36 degree
  Celsius but the temperature that is used is 34 degrees.
* To reproduce Figure 12C, the "presence of dendritic shunt conductances" meant
  setting g_L to 0.15 mS/cm^2 for the whole neuron.
* Other small discrepancies with the paper -- values from the NEURON code were
  used whenever different from the values stated in the paper

::

    
    from brian2 import *
    from brian2.units.constants import (zero_celsius, faraday_constant as F,
                                        gas_constant as R)
    
    defaultclock.dt = 0.01*ms
    
    VT = -52*mV
    El = -76.5*mV  # from code, text says: -69.85*mV
    
    E_Na = 50*mV
    E_K = -100*mV
    C_d = 7.954  # dendritic correction factor
    
    T = 34*kelvin + zero_celsius # 34 degC (current-clamp experiments)
    tadj_HH = 3.0**((34-36)/10.0)  # temperature adjustment for Na & K (original recordings at 36 degC)
    tadj_m_T = 2.5**((34-24)/10.0)
    tadj_h_T = 2.5**((34-24)/10.0)
    
    shift_I_T = -1*mV
    
    gamma = F/(R*T)  # R=gas constant, F=Faraday constant
    Z_Ca = 2  # Valence of Calcium ions
    Ca_i = 240*nM  # intracellular Calcium concentration
    Ca_o = 2*mM  # extracellular Calcium concentration
    
    eqs = Equations('''
    Im = gl*(El-v) - I_Na - I_K - I_T: amp/meter**2
    I_inj : amp (point current)
    gl : siemens/meter**2
    
    # HH-type currents for spike initiation
    g_Na : siemens/meter**2
    g_K : siemens/meter**2
    I_Na = g_Na * m**3 * h * (v-E_Na) : amp/meter**2
    I_K = g_K * n**4 * (v-E_K) : amp/meter**2
    v2 = v - VT : volt  # shifted membrane potential (Traub convention)
    dm/dt = (0.32*(mV**-1)*(13.*mV-v2)/
            (exp((13.*mV-v2)/(4.*mV))-1.)*(1-m)-0.28*(mV**-1)*(v2-40.*mV)/
            (exp((v2-40.*mV)/(5.*mV))-1.)*m) / ms * tadj_HH: 1
    dn/dt = (0.032*(mV**-1)*(15.*mV-v2)/
            (exp((15.*mV-v2)/(5.*mV))-1.)*(1.-n)-.5*exp((10.*mV-v2)/(40.*mV))*n) / ms * tadj_HH: 1
    dh/dt = (0.128*exp((17.*mV-v2)/(18.*mV))*(1.-h)-4./(1+exp((40.*mV-v2)/(5.*mV)))*h) / ms * tadj_HH: 1
    
    # Low-threshold Calcium current (I_T)  -- nonlinear function of voltage
    I_T = P_Ca * m_T**2*h_T * G_Ca : amp/meter**2
    P_Ca : meter/second  # maximum Permeability to Calcium
    G_Ca = Z_Ca**2*F*v*gamma*(Ca_i - Ca_o*exp(-Z_Ca*gamma*v))/(1 - exp(-Z_Ca*gamma*v)) : coulomb/meter**3
    dm_T/dt = -(m_T - m_T_inf)/tau_m_T : 1
    dh_T/dt = -(h_T - h_T_inf)/tau_h_T : 1
    m_T_inf = 1/(1 + exp(-(v/mV + 56)/6.2)) : 1
    h_T_inf = 1/(1 + exp((v/mV + 80)/4)) : 1
    tau_m_T = (0.612 + 1.0/(exp(-(v/mV + 131)/16.7) + exp((v/mV + 15.8)/18.2))) * ms / tadj_m_T: second
    tau_h_T = (int(v<-81*mV) * exp((v/mV + 466)/66.6) +
               int(v>=-81*mV) * (28 + exp(-(v/mV + 21)/10.5))) * ms / tadj_h_T: second
    ''')
    
    # Simplified three-compartment morphology
    morpho = Cylinder(x=[0, 38.42]*um, diameter=26*um)
    morpho.dend = Cylinder(x=[0, 12.49]*um, diameter=10.28*um)
    morpho.dend.distal = Cylinder(x=[0, 84.67]*um, diameter=8.5*um)
    neuron = SpatialNeuron(morpho, eqs, Cm=0.88*uF/cm**2, Ri=173*ohm*cm,
                           method='exponential_euler')
    
    neuron.v = -74*mV
    # Only the soma has Na/K channels
    neuron.main.g_Na = 100*msiemens/cm**2
    neuron.main.g_K = 100*msiemens/cm**2
    # Apply the correction factor to the dendrites
    
    neuron.dend.Cm *= C_d
    neuron.m_T = 'm_T_inf'
    neuron.h_T = 'h_T_inf'
    
    mon = StateMonitor(neuron, ['v'], record=True)
    
    store('initial state')
    
    
    def do_experiment(currents, somatic_density, dendritic_density,
                      dendritic_conductance=0.0379*msiemens/cm**2,
                      HH_currents=True):
        restore('initial state')
        voltages = []
        neuron.P_Ca = somatic_density
        neuron.dend.distal.P_Ca = dendritic_density * C_d
        # dendritic conductance (shunting conductance used for Fig 12C)
        neuron.gl = dendritic_conductance
        neuron.dend.gl = dendritic_conductance * C_d
        if not HH_currents:
            # Shut off spiking (for Figures 12B and 12C)
            neuron.g_Na = 0*msiemens/cm**2
            neuron.g_K = 0*msiemens/cm**2
        run(180*ms)
        store('before current')
        for current in currents:
            restore('before current')
            neuron.main.I_inj = current
            print('.', end='')
            run(320*ms)
            voltages.append(mon[morpho].v[:])  # somatic voltage
        return voltages
    
    
    ## Run the various variants of the model to reproduce Figure 12
    mpl.rcParams['lines.markersize'] = 3.0
    fig, axes = plt.subplots(2, 2)
    print('Running experiments for Figure A1 ', end='')
    voltages = do_experiment([50, 75]*pA, somatic_density=1.7e-5*cm/second,
                             dendritic_density=1.7e-5*cm/second)
    print(' done.')
    cut_off = 100*ms  # Do not display first part of simulation
    axes[0, 0].plot((mon.t - cut_off) / ms, voltages[0] / mV, color='gray')
    axes[0, 0].plot((mon.t - cut_off) / ms, voltages[1] / mV, color='black')
    axes[0, 0].set(xlim=(0, 400), ylim=(-80, 40), xticks=[],
                   title='A1: Uniform T-current density', ylabel='Voltage (mV)')
    axes[0, 0].spines['right'].set_visible(False)
    axes[0, 0].spines['top'].set_visible(False)
    axes[0, 0].spines['bottom'].set_visible(False)
    
    print('Running experiments for Figure A2 ', end='')
    voltages = do_experiment([50, 75]*pA, somatic_density=1.7e-5*cm/second,
                             dendritic_density=9.5e-5*cm/second)
    print(' done.')
    cut_off = 100*ms  # Do not display first part of simulation
    axes[1, 0].plot((mon.t - cut_off) / ms, voltages[0] / mV, color='gray')
    axes[1, 0].plot((mon.t - cut_off) / ms, voltages[1] / mV, color='black')
    axes[1, 0].set(xlim=(0, 400), ylim=(-80, 40),
                   title='A2: High T-current density in dendrites',
                   xlabel='Time (ms)', ylabel='Voltage (mV)')
    axes[1, 0].spines['right'].set_visible(False)
    axes[1, 0].spines['top'].set_visible(False)
    
    print('Running experiments for Figure B ', end='')
    currents = np.linspace(0, 200, 41)*pA
    voltages_somatic = do_experiment(currents, somatic_density=56.36e-5*cm/second,
                                     dendritic_density=0*cm/second,
                                     HH_currents=False)
    voltages_somatic_dendritic = do_experiment(currents, somatic_density=1.7e-5*cm/second,
                                               dendritic_density=9.5e-5*cm/second,
                                               HH_currents=False)
    print(' done.')
    maxima_somatic = Quantity(voltages_somatic).max(axis=1)
    maxima_somatic_dendritic = Quantity(voltages_somatic_dendritic).max(axis=1)
    axes[0, 1].yaxis.tick_right()
    axes[0, 1].plot(currents/pA, maxima_somatic/mV,
                    'o-', color='black', label='Somatic only')
    axes[0, 1].plot(currents/pA, maxima_somatic_dendritic/mV,
                    's-', color='black', label='Somatic & dendritic')
    axes[0, 1].set(xlabel='Injected current (pA)', ylabel='Peak LTS (mV)',
                   ylim=(-80, 0))
    axes[0, 1].legend(loc='best', frameon=False)
    
    print('Running experiments for Figure C ', end='')
    currents = np.linspace(200, 400, 41)*pA
    voltages_somatic = do_experiment(currents, somatic_density=56.36e-5*cm/second,
                                     dendritic_density=0*cm/second,
                                     dendritic_conductance=0.15*msiemens/cm**2,
                                     HH_currents=False)
    voltages_somatic_dendritic = do_experiment(currents, somatic_density=1.7e-5*cm/second,
                                               dendritic_density=9.5e-5*cm/second,
                                               dendritic_conductance=0.15*msiemens/cm**2,
                                               HH_currents=False)
    print(' done.')
    maxima_somatic = Quantity(voltages_somatic).max(axis=1)
    maxima_somatic_dendritic = Quantity(voltages_somatic_dendritic).max(axis=1)
    axes[1, 1].yaxis.tick_right()
    axes[1, 1].plot(currents/pA, maxima_somatic/mV,
                    'o-', color='black', label='Somatic only')
    axes[1, 1].plot(currents/pA, maxima_somatic_dendritic/mV,
                    's-', color='black', label='Somatic & dendritic')
    axes[1, 1].set(xlabel='Injected current (pA)', ylabel='Peak LTS (mV)',
                   ylim=(-80, 0))
    axes[1, 1].legend(loc='best', frameon=False)
    
    plt.tight_layout()
    plt.show()
    

