XPP version of RS cell. Best fit of Fig1c. Documentation.

NeuroML2/README.md

@@ -0,0 +1,7 @@
+This folder contains the NeuroML2 versions of the models. The folders contain the following:
+
+- cells: The 5 cells that were provided with the original NEURON code, implemented in NeuroML2
+- channels: All the ion channels used by the 5 cells
+- figures: Figures from the paper that have been reproduced
+
+To run any of these models, find the LEMS_XYZ.xml files and run them via `jnml LEMS_XYZ.xml`. See [Installation Instructions](https://github.com/OpenSourceBrain/PospischilEtAl2008) for more details.

README.md

@@ -1,6 +1,62 @@
-### Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons
 
-Conversion to NeuroML of cell models from: [Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons](http://link.springer.com/article/10.1007/s00422-008-0263-8), Martin Pospischil, Maria Toledo-Rodriguez, Cyril Monier, Zuzanna Piwkowska, Thierry Bal, Yves Frégnac, Henry Markram and Alain Destexhe, Biological Cybernetics, 2008.
+## Minimal Hodgkin–Huxley type models for four common classes of cortical and thalamic neurons
+
+[![Build Status](https://travis-ci.org/OpenSourceBrain/PospischilEtAl2008.svg?branch=master)](https://travis-ci.org/OpenSourceBrain/PospischilEtAl2008)
+
+#### Overview of the Model
+
+[Pospischil et. al. 2008](http://link.springer.com/article/10.1007/s00422-008-0263-8) describes conductance-based ([Hodgin-Huxley](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model)) models of four different classes of [cortical](https://en.wikipedia.org/wiki/Cerebral_cortex) and [thalamic](https://en.wikipedia.org/wiki/Thalamus) neurons. The classes represented are:
+
+- Regular Spiking (RS) cells
+- Fast Spiking (FS) cells
+- Intrinsicly Bursting (IB)
+- Low-Threshold Spike (LTS) cells
+
+All cell models are composed of [ion channel](https://en.wikipedia.org/wiki/Ion_channel) models generating the following currents:
+
+- I<sub>Na</sub> Voltage dependent Na<sup>+</sup> current for depolarization phase of action potentials (APs)
+- I<sub>Kd</sub> Delayed rectifier K<sup>+</sup> current for hyperpolarization phase of APs
+- I<sub>M</sub> Slow non-inactivating K<sup>+</sup> current for spike-frequency adaptation
+- I<sub>L</sub> High-threshold Ca<sup>2+</sup> current for burst generation
+- I<sub>T</sub> Low-threshold Ca<sup>2+</sup> current for rebound-burst generation
+
+Leak and input current models are present in all cells as well.
+
+These models demonstrate that a wide variety of spiking behaviors can be implemented using just a few types of ion channels. 
+
+**Original Reference:**
+
+[Minimal Hodgkin–Huxley type models for different classes of cortical and thalamic neurons](http://link.springer.com/article/10.1007/s00422-008-0263-8), Martin Pospischil, Maria Toledo-Rodriguez, Cyril Monier, Zuzanna Piwkowska, Thierry Bal, Yves Frégnac, Henry Markram and Alain Destexhe, *Biological Cybernetics*, 2008.
+
+#### Demonstrated Physiological Properties
+
+Specifically, these models demonstrate the following properties:
+
+**Regular Spiking**
+The use of I<sub>Na</sub> and I<sub>Kd</sub> currents and applying a constant input current is sufficient to generate a train of spikes with constant frequency. 
+
+![RS](Regular Spiking.png)
+
+**Frequency Adaptation**
+In addition to I<sub>Na</sub> and I<sub>Kd</sub> currents, including the I<sub>M</sub> current will result in a train of spikes where the inter-spike interval increases (and frequency decreases) with each spike in the train.
+
+![FA](Frequency adaptation.png)
+
+**Intrinsic Bursting**
+
+When a calcium-dependent I<sub>L</sub> current is added to a cell with I<sub>Na</sub>, I<sub>Kd</sub>, and I<sub>M</sub> currents, the cell membrane potential will display an initial high-frequency burst of spikes and then settle onto a regular, frequency-adapting spiking behavior. The I<sub>L</sub> current makes use of variable Ca<sup>2+</sup> reversal potential which is computed from the [Nernst equation](https://en.wikipedia.org/wiki/Nernst_equation) by tracking changes in the intracellular Ca<sup>2+</sup> concentration.
+
+![IB](Intrinsic bursting.png)
+
+**Rebound Bursting**
+
+Similarly, using the calcium-dependent I<sub>T</sub> current instead of I<sub>L</sub> results in a cell that will fire a series of high-frequency spikes after a negative (hyperpolarizing) current is withdrawn from cell. 
+
+![Rebound](Rebound.png)
+
+Figure produced with NEURON_ORIG LTS Cell with gcabar_it = 0.0012 and input current = -0.1 
+
+### Model Versions
 
 The original version of this model as uploaded to ModelDB can be found in directory [NEURON_ORIG](https://github.com/OpenSourceBrain/PospischilEtAl2008/tree/master/NEURON_ORIG). 
 An updated version of these scripts (in directory [NEURON_MODIFIED](https://github.com/OpenSourceBrain/PospischilEtAl2008/tree/master/NEURON_MODIFIED); 
@@ -9,8 +65,32 @@ An updated version of these scripts (in directory [NEURON_MODIFIED](https://gith
 A [neuroConstruct](http://www.neuroconstruct.org/) project [is included](https://github.com/OpenSourceBrain/PospischilEtAl2008/tree/master/neuroConstruct) which uses some of the ion channels to create cell models, but the 
 [NeuroML2 version](https://github.com/OpenSourceBrain/PospischilEtAl2008/tree/master/NeuroML2) of the model is much more complete.
 
-[![Build Status](https://travis-ci.org/OpenSourceBrain/PospischilEtAl2008.svg?branch=master)](https://travis-ci.org/OpenSourceBrain/PospischilEtAl2008)
+#### Why Convert to NeuroML?
+
+The original models were implemented for the [NEURON simulator](https://www.neuron.yale.edu/neuron/). We have converted the model to NeuroML. The modular, XML nature of NeuroML allows to quickly re-use this model in network simulations and our tools allow [automated conversion to other supported simulator formats](https://neuroml.org/mappings).
+
+#### NeuroML Version
+
+All five cells have been implemented in NeuroML and match the output produced by the cells provided in the [original NEURON code](http://senselab.med.yale.edu/ModelDB/ShowModel.cshtml?model=123623).
+
+#### XPP Version
+
+The Regular Spiking (RS) cell was also implemented in [XPP simulator](http://www.math.pitt.edu/~bard/xpp/xpp.html) to enable rapid debugging of parameter and other issues with the original models.
+
+
+### Installation Instructions
+
+1. [Download the Model Files](archive/master.zip), or clone the repository using git: `git clone https://github.com/OpenSourceBrain/PospischilEtAl2008.git`
+2. [Follow instructions to Install jNeuroML](https://github.com/NeuroML/jNeuroML) for the **jnml** executable. On Windows, you may also need [SVN](https://subversion.apache.org/packages.html#windows). Alternatively install [PyNeuroML](https://github.com/NeuroML/pyNeuroML) for the **pynml** executable. 
+3. Set the $PATH and $JNML_HOME variables as described in [#2](https://github.com/NeuroML/jNeuroML)
+4. Extract the model files to a folder.
+5. For Figure 1: Change to NeuroML2/Figure1 folder. Type `jnml LEMS_Figure1C.xml` or `pynml LEMS_Figure1C.xml`.
+6. For RS (and other cells): Change to NeuroML2/cells/RS (or other) folder. Type `jnml LEMS_RS.xml` or `pynml LEMS_RS.xml`. 
+7. Windows with the plotted figures should show up as can be seen above.
+
 
+### Issues
 
+See [Issues Tab](issues) for reported issues with the original model and the conversion. 
 
 

XPP/README.md

@@ -0,0 +1,5 @@
+This folder contains the XPP version of the RS (Regular Spiking) cell implementing Figure 1C from the paper. There are files for:
+
+ - _fromFigure: This is the model that is using the parameters for Figure 1C as **described in the paper figure**
+ - _fromNEURON: Parameters as **described in the original NEURON code**
+ - _bestFit: A combination of parameters from paper, NEURON, **and as discovered through parameter search to find the best fit to the published figure**

XPP/RS_bestFit.ode

@@ -0,0 +1,118 @@
+# Pospischil, et. al. Regular Spiking (RS) cell model
+# Figure 1c
+# Parameters that have been found to best reproduce the published figure
+# *** indicates discovered parameter
+
+# ---- PARAMETERS ---- #
+
+# 1 nF in section 2.1 & NEURON
+par Cm=1
+
+# Both 96 um in figure & NEURON
+par Ldiam = 0.96
+
+# from figure & NEURON
+par Eleak=-70
+
+# from section 2.2.1 Na current & NEURON
+par Ena=50
+
+# -90mV in section 2.2.1 K current
+# -100mV in NEURON hoc file
+# -100mv best fit to figure
+par Ek=-100
+
+# ??? Not in paper or figure
+# -63mV in NEURON mod file
+# -55mV in NEURON hoc file
+# -68.6mV best fit to figure ***
+par Vt=-68.6
+
+# 1e-4  S/cm2 in figure & NEURON
+par gleak=0.1
+
+# 0.05  S/cm2 in figure & NEURON
+par gna=50
+
+# 0.005 S/cm2 in figure & NEURON
+par gkd=5
+
+# 7e-5  S/cm2 in figure & NEURON
+par gm=0.07
+
+# 4000ms in section 2.2.2
+# 1000ms in NEURON hoc file
+# 1000ms best fit to figure
+par tmax=1000
+
+# Estimated by pixel count del=200, dur=500
+# In NEURON hoc del=300, dur=400
+par Idelay=200
+par Iduration=500
+
+# 0.5nA in figure
+# 0.75nA in NEURON hoc
+# 1.87 best fit to figure ***
+par I=1.87
+
+# No DC current in paper or NEURON
+# Figure starts and stays at -65mV
+# DC current of 1.12 needed to keep at -65mV ***
+par Idc=1.12
+
+# initial conditions
+# Figure starts at -65mV
+# -70mV in NEURON hoc
+init V=-65
+
+
+# ---- EQUATIONS ---- #
+
+# Lateral area
+area=2*PI*Ldiam/2*Ldiam
+
+# Total Capacitance
+# 0.29 nF in figure
+# 0.28953 nF in NEURON
+Cap=area*Cm
+
+# Current balance
+V'=(-Ileak-INa-IKd-Im+Iinput/area)/Cap
+
+# Input current with DC, delay, and duration
+Iinput=if(t>=Idelay)then(if(t<=(Idelay+Iduration))then(I)else(0))else(0)+Idc
+
+# Leak current
+Ileak=gleak*(V-Eleak)
+
+# Na current
+INa=gna*m*m*m*h*(V-Ena)
+
+m'=am*(1-m)-bm*m
+h'=ah*(1-h)-bh*h
+
+am=(-0.32*(V-Vt-13))/(exp(-(V-Vt-13)/4)-1)
+bm=(0.28*(V-Vt-40))/(exp((V-Vt-40)/5)-1)
+ah=0.128*exp(-(V-Vt-17)/18)
+bh=(4)/(1+exp(-(V-Vt-40)/5))
+
+# Kd current
+IKd=gkd*n*n*n*n*(V-Ek)
+
+n'=an*(1-n)-bn*n
+
+an=(-0.032*(V-Vt-15))/(exp(-(V-Vt-15)/5)-1)
+bn=0.5*exp(-(V-Vt-10)/40)
+
+# Km current
+Im=gm*p*(V-Ek)
+
+p'=(pinf-p)/tp
+
+pinf=1/(1+exp(-(V+35)/10))
+tp=tmax/(3.3*exp((V+35)/20)+exp(-(V+35)/20))
+
+# integrator params
+@ maxstor=80000,total=1000,bound=10000,xlo=0,xhi=1000,ylo=-75,yhi=40
+@ meth=cvode,atol=0.0001,toler=0.0001,dt=0.3
+done

XPP/RS_fromFigure.ode

@@ -0,0 +1,112 @@
+# Pospischil, et. al. Regular Spiking (RS) cell model
+# Figure 1c
+
+# ---- PARAMETERS ---- #
+
+# 1 nF in section 2.1 & NEURON
+par Cm=1
+
+# Both 96 um in figure & NEURON
+par Ldiam = 0.96
+
+# from figure & NEURON
+par Eleak=-70
+
+# from section 2.2.1 Na current & NEURON
+par Ena=50
+
+# -90mV in section 2.2.1 K current
+# -100mV in NEURON hoc file
+par Ek=-90
+
+# ??? Not in paper or figure
+# -63mV in NEURON mod file
+# -55mV in NEURON hoc file
+par Vt=-63
+
+# 1e-4  S/cm2 in figure & NEURON
+par gleak=0.1
+
+# 0.05  S/cm2 in figure & NEURON
+par gna=50
+
+# 0.005 S/cm2 in figure & NEURON
+par gkd=5
+
+# 7e-5  S/cm2 in figure & NEURON
+par gm=0.07
+
+# 4000ms in section 2.2.2
+# 1000ms in NEURON hoc file
+par tmax=4000
+
+# Estimated by pixel count del=200, dur=500
+# In NEURON hoc del=300, dur=400
+par Idelay=200
+par Iduration=500
+
+# 0.5nA in figure
+# 0.75nA in NEURON hoc
+par I=5
+
+# No DC current in paper or NEURON
+# Figure starts and stays at -65mV
+# DC current of 0.516 needed to keep at -65mV
+par Idc=0
+
+# initial conditions
+# Figure starts at -65mV
+# -70mV in NEURON hoc
+init V=-65
+
+
+# ---- EQUATIONS ---- #
+
+# Lateral area
+area=2*PI*Ldiam/2*Ldiam
+
+# Total Capacitance
+# 0.29 nF in figure
+# 0.28953 nF in NEURON
+Cap=area*Cm
+
+# Current balance
+V'=(-Ileak-INa-IKd-Im+Iinput/area)/Cap
+
+# Input current with DC, delay, and duration
+Iinput=if(t>=Idelay)then(if(t<=(Idelay+Iduration))then(I)else(0))else(0)+Idc
+
+# Leak current
+Ileak=gleak*(V-Eleak)
+
+# Na current
+INa=gna*m*m*m*h*(V-Ena)
+
+m'=am*(1-m)-bm*m
+h'=ah*(1-h)-bh*h
+
+am=(-0.32*(V-Vt-13))/(exp(-(V-Vt-13)/4)-1)
+bm=(0.28*(V-Vt-40))/(exp((V-Vt-40)/5)-1)
+ah=0.128*exp(-(V-Vt-17)/18)
+bh=(4)/(1+exp(-(V-Vt-40)/5))
+
+# Kd current
+IKd=gkd*n*n*n*n*(V-Ek)
+
+n'=an*(1-n)-bn*n
+
+an=(-0.032*(V-Vt-15))/(exp(-(V-Vt-15)/5)-1)
+bn=0.5*exp(-(V-Vt-10)/40)
+
+# Km current
+Im=gm*p*(V-Ek)
+
+p'=(pinf-p)/tp
+
+pinf=1/(1+exp(-(V+35)/10))
+tp=tmax/(3.3*exp((V+35)/20)+exp(-(V+35)/20))
+
+# integrator params
+@ maxstor=80000,total=1000,bound=10000,xlo=0,xhi=1000,ylo=-75,yhi=40
+@ meth=cvode,atol=0.0001,toler=0.0001,dt=0.3
+done