Optical stimulation of a channelrhodopsin-2 positive pyramidal neuron model (Foutz et al 2012)

 Download zip file   Auto-launch 
Help downloading and running models
A computational tool to explore the underlying principles of optogenetic neural stimulation. This "light-neuron" model consists of theoretical representations of the light dynamics generated by a fiber optic in brain tissue, coupled to a multicompartment cable model of a cortical pyramidal neuron (Hu et al. 2009, ModelDB #123897) embedded with channelrhodopsin-2 (ChR2) membrane dynamics. Simulations predict that the activation threshold is sensitive to many of the properties of ChR2 (density, conductivity, and kinetics), tissue medium (scattering and absorbance), and the fiber-optic light source (diameter and numerical aperture). This model system represents a scientific instrument to characterize the effects of optogenetic neuromodulation, as well as an engineering design tool to help guide future development of optogenetic technology.
1 . Foutz TJ, Arlow RL, McIntyre CC (2012) Theoretical principles underlying optical stimulation of a channelrhodopsin-2 positive pyramidal neuron. J Neurophysiol 107:3235-45 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell; Axon; Channel/Receptor; Dendrite;
Brain Region(s)/Organism:
Cell Type(s): Neocortex L5/6 pyramidal GLU cell;
Channel(s): I K; I M; I K,Ca; I Sodium; I Calcium; Channelrhodopsin (ChR);
Gap Junctions:
Simulation Environment: NEURON; Python;
Model Concept(s): Detailed Neuronal Models; Action Potentials; Parameter sensitivity; Intracortical Microstimulation; Neuromodulation;
Implementer(s): MacIntyre, CC [ccm4 at case.edu]; Foutz, Thomas J [tfoutz at uw.edu]; Arlow, Richard L [richard.arlow at case.edu];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; I K; I M; I K,Ca; I Sodium; I Calcium; Channelrhodopsin (ChR);

26 Ago 2002 Modification of original channel to allow variable time
step and to correct an initialization error.
    Done by Michael Hines(michael.hines@yale.e) and Ruggero
Scorcioni(rscorcio@gmu.edu) at EU Advance Course in Computational
Neuroscience. Obidos, Portugal
11 Jan 2007
    Glitch in trap where (v/th) was where (v-th)/q is. (thanks Ronald
van Elburg!)


Sodium channel, Hodgkin-Huxley style kinetics.

Kinetics were fit to data from Huguenard et al. (1988) and Hamill et
al. (1991)

qi is not well constrained by the data, since there are no points
between -80 and -55.  So this was fixed at 5 while the thi1,thi2,Rg,Rd
were optimized using a simplex least square proc

voltage dependencies are shifted approximately from the best
fit to give higher threshold

Author: Zach Mainen, Salk Institute, 1994, zach@salk.edu

20150512 NTC
Cleaned code up a bit.
Made threadsafe and replaced vtrap0 with efun,
which is a better approximation in the vicinity of a singularity.
See na.mod in ModelDB entry 2488.

Special comment:

This mechanism borrowed heavily from na.mod in ModelDB entry 2488.
That code, which was intended to be used only at 37 deg C 
(see comments from 20120514 in that file),
calculated ionic conductance as the product
  g = tadj*gbar*product_of_gating_variables
  tadj = q10^((celsius - temp)/10)
  temp is the "reference temperature" (at which the gating variable
    time constants were originally determined)
  celsius is the "operating temperature"
This deviates from the standard HH formula
  g = gbar*product_of_gating_variables
and has the unfortunate consequence of not only making the 
effective channel density differ from the nominal (i.e. user-assigned) 
channel density, but it would also make the effective channel density
depend on temperature!
Sooner or later this is guaranteed to confound studies of the effects 
of temperature on model operation.
It would also be a debugging nightmare, not least because
the ModelView tool--so handy for discovering the properties of a 
model and verifying a close match between the computational model
and the conceptual model--would report the NOMINAL channel density,
not the effective channel density.

To eliminate this problem, every statement of the form
  g = tadj*gbar*gating variables
has been replaced by 
  g = gbar*gating variables
Furthermore, the numerical value assigned to gbar,
whether by an assignment statement in the PARAMETER block
or by a hoc or Python statement, 
must now use the actual conductance density.
For this particular mechanism, tadj at 37 deg C is 3.20936
so the actual conductance density is 3.20936 times 
the nominal conductance density.

  USEION na READ ena WRITE ina
  RANGE m, h, gna, gbar
  GLOBAL tha, thi1, thi2, qa, qi, qinf, thinf
  RANGE minf, hinf, mtau, htau
  GLOBAL Ra, Rb, Rd, Rg
  GLOBAL q10, temp, tadj, vmin, vmax, vshift

  (mA) = (milliamp)
  (mV) = (millivolt)
  (pS) = (picosiemens)
  (um) = (micron)

: this assignment doesn't matter
: if it is overridden by assignments in hoc or Python
  gbar = 1000 (pS/um2) : 0.1 mho/cm2
  vshift = -5 (mV)    : voltage shift (affects all)

  tha  = -43  (mV)    : v 1/2 for act    (-42)
  qa   = 7  (mV)      : act slope
  Ra   = 0.182  (/ms) : open (v)
  Rb   = 0.124  (/ms) : close (v)

  thi1  = -50  (mV)   : v 1/2 for inact
  thi2  = -75  (mV)   : v 1/2 for inact
  qi   = 5  (mV)      : inact tau slope
  thinf  = -72  (mV)  : inact inf slope
  qinf  = 6.2  (mV)   : inact inf slope
  Rg   = 0.0091  (/ms) : inact (v)
  Rd   = 0.024  (/ms)  : inact recov (v)

  temp = 23  (degC)   : original temp
  q10  = 2.3          : temperature sensitivity

:  v     (mV)
:  dt    (ms)
:  celsius    (degC)
  vmin = -120  (mV)
  vmax = 100  (mV)

  v      (mV)
  celsius (degC)
  a      (/ms)
  b      (/ms)
  ina    (mA/cm2)
  gna    (pS/um2)
  ena    (mV)
  minf   hinf
  mtau (ms)  htau (ms)

STATE { m h }

  : since tadj is a per-thread GLOBAL
  : all threads must calculate its value at initialization
  tadj = q10^((celsius - temp)/(10 (degC)))
  m = minf
  h = hinf

  SOLVE states METHOD cnexp
:  gna = tadj*gbar*m*m*m*h
  gna = gbar*m*m*m*h
  ina = (1e-4) * gna * (v - ena)

: LOCAL mexp, hexp

DERIVATIVE states {   :Computes state variables m and h
  trates(v+vshift)    :  at the current v
  m' =  (minf-m)/mtau
  h' =  (hinf-h)/htau
:  m = m + mexp*(minf-m)
:  h = h + hexp*(hinf-h)

PROCEDURE trates(v (mV)) {
  TABLE minf, hinf, mtau, htau
  DEPEND celsius, temp, Ra, Rb, Rd, Rg, tha, thi1, thi2, qa, qi, qinf
  FROM vmin TO vmax WITH 199

  rates(v): not consistently executed from here if usetable == 1

:  tadj = q10^((celsius - temp)/10)
:  tinc = -dt * tadj

:  mexp = 1 - exp(tinc/mtau)
:  hexp = 1 - exp(tinc/htau)

PROCEDURE rates(vm (mV)) {
:  LOCAL  a, b

:  a = trap0(vm,tha,Ra,qa)
:  b = trap0(-vm,-tha,Rb,qa)
  a = Ra*qa*efun((tha-vm)/qa)
  b = Rb*qa*efun((vm-tha)/qa)

  tadj = q10^((celsius - temp)/(10 (degC)))

  mtau = 1/tadj/(a+b)
  minf = a/(a+b)

:  mtau = 1/(a+b)
:  minf = a*mtau

  :"h" inactivation

:  a = trap0(vm,thi1,Rd,qi)
:  b = trap0(-vm,-thi2,Rg,qi)
  a = Rd*qi*efun((thi1-vm)/qi)
  b = Rg*qi*efun((vm-thi2)/qi)

  htau = 1/tadj/(a+b)
  hinf = 1/(1+exp((vm-thinf)/qinf))

:  htau = 1/(a+b)
:  hinf = 1/(1+exp((vm-thinf)/qinf))

FUNCTION trap0(v,th,a,q) {
:  if (fabs(v/th) > 1e-6) {
  if (fabs((v-th)/q) > 1e-6) {
    trap0 = a * (v - th) / (1 - exp(-(v - th)/q))
  } else {
    trap0 = a * q

FUNCTION efun(z) {
  if (fabs(z) < 1e-6) {
    efun = 1 - z/2
    efun = z/(exp(z) - 1)