Known issues with Traub et al 2005.¶
This is a quite complex and detailed model and as discussed in the original paper
Any model, even of a small bit of cortex, is subject to difficulties and hazards: limited data, large numbers of parameters, criticisms that models with complexity comparable to the modeled system cannot be scientifically useful, the expense and slowness of the necessary computations, and serious uncertainties as to how a complex model can be compared with experiment and shown to be predictive.
The above difficulties and hazards are too real to be dismissed readily. In our opinion, the only way to proceed is through a state of denial that any of the difficulties need be fatal. The reader must then judge whether the results, preliminary as they must be, help our understanding.
Even the published Fortran version of this model was acknowledged to be incomplete. Each conversion of this model will deviate to a small or large extent from this version.
Questions about physiological properties of model¶
Dependence on Fast Regular Bursting cells for oscillatory behaviour
Prevalence of gap junctions
High current threshold for deep pyramidal firing
Not tested with external synaptic input
Limitations of the conversion of the model to NEURON¶
It is useful to read the notes on conversion of this model to NEURON from Fortran by Tom Morse and Michael Hines
Slightly different method of running the simulation (e.g. in Neuron information about spike is sent immediately, in Fortran every 0.1 ms )
**Sum of transmembrane currents in every single cells sums up to 0 only if you use cvode_active*
**Diffrent behaviour of NMDA synapse when thalamus is disconnected* (some bug in Fortran version?)
In Fortran code:
z = 0.d0 ! thalamus disconnected gAMPA_TCR_to_suppyrRS = z * gAMPA_TCR_to_suppyrRS gNMDA_TCR_to_suppyrRS = z * gNMDA_TCR_to_suppyrRS gAMPA_TCR_to_suppyrFRB = z * gAMPA_TCR_to_suppyrFRB gNMDA_TCR_to_suppyrFRB = z * gNMDA_TCR_to_suppyrFRB ...
gNMDA_TCR_to_suppyrFRB becomes 0. Then when you compute NMDA activation
from TCR to suppyrFRB
.... ! NMDA part if (delta.le.5.d0) then gNMDA_suppyrFRB(k,L) = gNMDA_suppyrFRB(k,L) + & gNMDA_TCR_to_suppyrFRB * delta * 0.2d0 else dexparg = (delta - 5.d0)/tauNMDA_TCR_to_suppyrFRB if (dexparg.le.5.d0) then z = dexptablesmall (int(dexparg*1000.d0)) else if (dexparg.le.100.d0) then z = dexptablebig (int(dexparg*10.d0)) else z = 0.d0 endif gNMDA_suppyrFRB(k,L) = gNMDA_suppyrFRB(k,L) + & gNMDA_TCR_to_suppyrFRB * z endif c Test for NMDA saturation z = NMDA_saturation_fact * gNMDA_TCR_to_suppyrFRB if (gNMDA_suppyrFRB(k,L).gt.z) & gNMDA_suppyrFRB(k,L) = z ! end NMDA part ....
It seems that this piece of code, more precisely the last three lines:
c Test for NMDA saturation z = NMDA_saturation_fact * gNMDA_TCR_to_suppyrFRB if (gNMDA_suppyrFRB(k,L).gt.z) & gNMDA_suppyrFRB(k,L) = z </pre> kills completely NMDA activation of suppyrFRB cells from all the other populations, not just TCR (except from nontuftRS cells, nontuftRS - suppyrFRB NMDA conductance is calculated after this block). In Neuron version there is no such behaviour. An *updated version* of this model in NEURON is being worked on "here":https://github.com/hglabska/Thalamocortical/tree/Neuron_version_simplified_groucho_file/Neuron. The version allows to modify easily the network, e.g. to add new population (version commited on 26 June 2013 and later), replace one template by another e.g. tuftIB Traub cell by "Hay cell":http://senselab.med.yale.edu/ModelDb/ShowModel.asp?model=139653 ( version commited on 04 July 2013 or later). The main groucho.hoc file is simpler and much shorter (about 10 times), parameters like AMPA, GABA, NMDA conductances, connections between populations are defined in separated files. h4. Tests for "Neuron":http://senselab.med.yale.edu/ModelDb/ShowModel.asp?model=82894 and "Fortran":https://github.com/hglabska/Thalamocortical/tree/master/Fortran_ifc version . Trying to reproduce results from the "article":http://www.ncbi.nlm.nih.gov/pubmed/15525801 *Remark 1* In Fortran version, compilation flag -finit-local-zero , seems to be important! *Remark 2* If you want to run the Fortran version locally on less than 14 cores you can do this in this way: <pre> echo localhost >> my_hostfile mpirun -np 14 --hostfile my_hostfile ./groucho
Thanks to kindness of Roger Traub, who sent us parameters which were used to generate figures 2. and 7. in the article , we were able to test how well we can reproduce the results on different version of the model.
Results from Appendix A - activity of single cells after applying some current to the soma, were reproduce reasonable well in Neuron version. For more data look here .
To compare the single cell result in Neuron with Fortran version you can use this code with makefile.single_cell instead of makefile. This version contains additional 14 programs to simulate single cell from every of 14 populations.
The biggest challenge in Appendix A is to reproduce fig A4C: applying some pulse current in apical dendrite caused somatic burst.
First difficulties is to estimate the amplitude of the current (It is not describe in article).
I =3* /10)) * /20)) nA,
looks reasonable well:
but Neuron result doesn’t look similar like the result in the article (colours: green D1, black D2, red soma):
also Fortran version (tuftIB.f ) of the model failed to reproduce the somatic burst with the same stimulus.
It is possible to obtain this somatic spikes ( in both Neuron and Fortran version) after depolarizing the soma by 1nA current and increasing the apical stimulus 3 times. Decreasing the depolarizing somatic currents two times (0.5 nA) , or using the apical stimulus like at the beginning (I =3* /10)) * /20)) ), caused that the somatic spikes disappear.
Remark : Look at the difference in the somatic membrane potential after the burst, between Fortran and Neuron versions.
For more data look here ( EPSP means apical stimulus with amplitude I =3* /10)) * /20))):
EPSP, 3 * EPSP
“Simulation of kainate-induced gamma oscillations”
The results in both Neuron and Fortran version looks quite similar. Only be aware that activity of suppyrRS differs much between individual cells. One questionable issue is appearance of the burst after about 1500 ms in Fortran and nearly 1200 ms in Neuron version (not shown here), which they didn’t report in the article.
You can download the data (+ rasterplot) for Fig 2 from Fortran and Neuron simulation: Fortran data and Neuron data.
For Neuron simulation you can also compare the result with simulation using the “traub_exact()” algorithm: Neuron traub_excat() data. More about “traub_excat()” algorithm you can read in notes on conversion of this model to NEURON from Fortran by Tom Morse and Michael Hines.
“Effects of disinhibition in model (cortex only, with thalamus disconnected), when there are open gap junctions between the axons of the respective principal cell populations (superficial pyramids, spiny stellates, layer 5 pyramids, layer 6 pyramids), and spiny stellates are strongly interconnected by AMPA receptors .”
Figure 7 from the article
In the article they raported about consisting of 17 burst complexes that terminate spontaneously. The last 5 of the bursts are shown. Results from the Fortran version are very similar but only 14 bursts appears. In Neuron version the result is much different.
You can download the data (+ rasterplot) for Fig 7A from Fortran and Neuron simulation: Fortran data and Neuron data.
For Neuron simulation you can also compare the result with simulation using the “traub_exact()” algoritm: Neuron traub_excat() data. More about “traub_excat()” algoritm you can read in notes on conversion of this model to NEURON from Fortran by Tom Morse and Michael Hines.
Response to simple stimulus - comparison between Fortran and Neuron versions.¶
Gap junctions are closed, thalamus is connected with cortex.
Stimulus: current injection to thalamic (TCR) somas . Current delay 300 ms, duration 2 ms, amplitude 3 nA.
In normal case there is small, short response in layers 2/3, 4 and inhibitory neurons in layers 5/6. The answer is much better visible if we decrease GABA conductances, but still there is no response in layer 5 and 6 in pyramidal cells (except ectopic spikes).
No response in the cortex in normal case. Answer in layers 2/3, 4 and inhibitory neurons in layers 5/6 after decreasing GABA conductances, but activity in layers 2/3 is shorter than in Fortran case, single spike in pyramidal cells layer 6 and no response in layer 5 (only ectopic spikes).
Applying additional current tu somas in pyramidal cells in layer 5 (1 nA) and 6 (0.75 nA) ( awake = 1 in Neuron version).
The additional stimulus is to big, a lot of spontaneous burst in every case. In Fortran version response in layer 2/3 lasts longer.
Additional current to somas; 0.5 nA in pyramidal cells in layer 5 and 0.375 nA in somas of pyramids in layer 6.
The additional stimulus is still to big in Fortran version ( a lot of spontaneous burst).
Spontaneous burst still exist but there are very seldom (not shown on the picture). Now can observe response in every layer.
Additional current to somas; 0.2 nA in pyramidal cells in layer 5 and 0.15 nA in somas of pyramids in layer 6.
All layers answer to stimulus, only response in pyramids layer 5 and 6 is quite late.
Again no response in layers 5 and single spike or no response in layer 6 pyramids.
- Fortran and Neuron code doesn’t generate the soma output even when gap junctions are closed
- in Fortran version response in layer 2/3 is more complex, 3 bursts versus 1 (why? )
- when gap junctions are closed in normal condition when GABA conductance is not decreased, response in layer 2/3 pyramids last extremely short (single spikes)
Limitations of the conversion of the model to MOOSE¶
Limitations of the conversion of the model to NeuroML¶
Optimal spatial discretisation for each cell needs to be investigated
Important details of the process of conversion of the cell models to NeuroML, and matching cell behaviour across simulators is present in the 2010 NeuroML paper.
The spatial discretisation of the cells influenced precise spike timing. Changing the number of compartments/points used to calculate the membrane potential changed the timing of the cell (e.g. changing the value of nseg in NEURON on all sections). See below for an example of how 3 cells with differing numbers of compartments converged at different rates. A) Nucleus reticularis thalami (nRT) cell; B) Superficial Low Threshold spiking (LTS) cell; C) Layer 6 Non-tufted Regular Spiking pyramidal cell. Traces for NEURON (black) and MOOSE (green) and GENESIS (red).
NMDA conductance wave form
The NMDA synapse model used in the network has an unconventional form, with a scaling factor rising lineally between 0 and 5ms, and decaying exponentially. This can probably be approximated by a double exponential synapse (coupled with v & [Mg] dependent blocking mechanism).
Firing rate vs. injected current of cells
Many of the cells show unusual F/I curves.
Support in NeuroML
All model elements from the neuroConstruct generated network can be exported to valid NeuroML v1.8.1.