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

Note: data generated by Helena Głąbska and Chaitanya Chintaluri at the Nencki Institute in Warsaw from the Neuron version of this model (used in the comments below) can be found here.

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

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.


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.

Single Cell

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):

Neuron

also Fortran version (tuftIB.f ) of the model failed to reproduce the somatic burst with the same stimulus.

Fortran

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.

Neuron

Fortran

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))*):

Fortran
EPSP

somatic current 1nA + 3*EPSP

somatic current 1nA + EPSP

somatic current 0.5 nA + 3*EPSP

Neuron
EPSP, 3 * EPSP

somatic current 1 nA + EPSP, 3*EPSP

somatic current 0.5 nA + EPSP, 3*EPSP

Figure 2

“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.

Figure 7

“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

7A

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.

7B
Results from the Fortran version looks again very similar, although gives much more complex bursts, at least 6, when prolong the simulation up to 2000 ms (results not shown here - download ) .

You can download the data (+ rasterplot) for Fig 7B from Fortran and Neuron simulation: Fortran data and Neuron data compare with Neuron traub_excat() data .

7C

You can download the data (+ rasterplot) for Fig 7C from Fortran and Neuron simulation. Fortran data and Neuron data compare with Neuron traub_excat() data .

7D

You can download the data (+ rasterplot) for Fig 7D from Fortran and Neuron simulation. Fortran data and Neuron data compare with Neuron traub_excat() data

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.

Fortran
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).

Neuron
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).

data

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.

data + code

Additional current to somas; 0.5 nA in pyramidal cells in layer 5 and 0.375 nA in somas of pyramids in layer 6.

Fortran
The additional stimulus is still to big in Fortran version ( a lot of spontaneous burst).

Neuron
Spontaneous burst still exist but there are very seldom (not shown on the picture). Now can observe response in every layer.

data

Additional current to somas; 0.2 nA in pyramidal cells in layer 5 and 0.15 nA in somas of pyramids in layer 6.

Fortran
All layers answer to stimulus, only response in pyramids layer 5 and 6 is quite late.

Neuron
Again no response in layers 5 and single spike or no response in layer 6 pyramids.

data

Conclusions/Remarks:

  • 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

TODO…

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.

Model can be exported to (mostly valid) NeuroML 2, but there is not yet an application that can handle such detailed NML2 models (but we’re working on it).

7paper.png (285 KB) Helena Głąbska, 19 Nov 2013 15:24

7A_small_labels.png (295 KB) Helena Głąbska, 19 Nov 2013 16:11

7B_small_labels.png (201 KB) Helena Głąbska, 21 Nov 2013 09:40

7C_small_labels.png (162 KB) Helena Głąbska, 21 Nov 2013 10:11

7D_small_labels.png (172 KB) Helena Głąbska, 21 Nov 2013 10:28

test2_labels.png (377 KB) Helena Głąbska, 22 Nov 2013 13:41

A4C.png (65.5 KB) Helena Głąbska, 25 Nov 2013 16:19

pulse.png (2.83 KB) Helena Głąbska, 25 Nov 2013 16:26

tuftIB_Neuron_voltage.png (3.03 KB) Helena Głąbska, 25 Nov 2013 16:29

voltage_tutfIB_fortran.png (14.4 KB) Helena Głąbska, 26 Nov 2013 12:45

fortran_burst.png (16.4 KB) Helena Głąbska, 26 Nov 2013 12:53

neuron_burst.png (3.72 KB) Helena Głąbska, 26 Nov 2013 13:30

thalamus_awake0.png (207 KB) Helena Głąbska, 05 Dec 2013 11:55

thalamus_awake1.png (236 KB) Helena Głąbska, 05 Dec 2013 11:55

thalamus_awake02.png (217 KB) Helena Głąbska, 05 Dec 2013 11:56

thalamus_awake05.png (222 KB) Helena Głąbska, 05 Dec 2013 11:56