Example: Destexhe_et_al_1998
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):
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()
