Single node with McIntyre equations fires without stimulus

[mzenker]

Hi,

I am new to NEURON and want to reproduce a model found in a paper by Mercadal et al (2017) http://stacks.iop.org/0031-9155/62/i=20/a=8060. They say they use a model for a myelinated axon where the behaviour of the nodes is described by the equations by McIntyre (2002). I have found the corresponding NEURON model in the ModelDB (No. 3810). In the paper, they say they just use the modified HH equations together with the cable equation by McNeal (1976). So I tried to make just a node with those equations, the next step would be to make more nodes and add an extracellular stimulus.

The problem I have is that the node is firing continuously without a stimulus, although I have initialized the potential to the value given in the model (-80 mV).

Here is the NEURON code I use. The parameters for the node are taken from the McIntyre model file.

Code: Select all

proc globals(){
	celsius=37			
	v_init=-80 //mV//  		
	dt=0.005 //ms//         	
	tstop=20
	nodeD=3.3
	nodelength=1.0
	rhoa=0.7e6 //Ohm-um//
}
globals()

create node
access node

proc initialize(){
	node {
		nseg=1
		diam=nodeD
		L=nodelength
		Ra=rhoa/10000
		cm=2
		insert axnode			
		//insert hh
	}
}
initialize()

axnode is the .mod file from the McIntyre model containing the modified HH equations. I just compiled it without change.
If I comment out the celsius statement, change v_init to -65 mV and use hh instead of axnode, I get the expected flat line at -65 mV.

[EDIT] I have added a stimulus (using the code from the McIntyre model) to verify that with hh, it fires if I give some current stimulus, and it doesn't if I set the current to zero. With axnode, it fires all the time, the shape of the first spike just changes a bit if I add a stimulus.[/EDIT]

What do I do wrong?

Thanks for a pointer or hint,

Matthias

[ted]

First, you need to be sure that your model's parameters and initialization are identical to those used by McIntyre et al.. Easiest way to check your model is to initialize your model, then use NEURON's Model View tool to examine parameter values. Easiest way to check the model used by McIntyre et al. is to get their model's source code from ModelDB
https://senselab.med.yale.edu/ModelDB/s ... model=3810
then initialize it and finally use NEURON's Model View tool to examine its parameter values.

[mzenker]

Hi Ted,

thank you for your answer.

I have now checked again that the Parameters are the same - as far as I use the McIntyre model, i.e. as far as the node and the axnode mechanism is concerned. Of course the model by McIntyre et al. is more complex, as it contains many nodes and also explicitly models the myelinated parts in between. But in the paper I try to reproduce, the authors claim that they don't do that. Here is how they describe their model:

The response of a nerve fiber to the voltage applied across the electrodes was determined using the cable model for a myelinated axon (McNeal 1976). Following this approach the TMV relative to the resting voltage (𝑉_𝑛=𝑉_𝑖,𝑛− 𝑉_𝑒,𝑛−𝑉_𝑟) at the nth node of Ranvier can be calculated by solving the following equation:
𝑑𝑉_𝑛/𝑑𝑡=1/𝐶_𝑚(𝐺𝑎(𝑉_(𝑛−1)−2𝑉_𝑛+𝑉_(𝑛+1)+𝑉_𝑒,(𝑛−1)−2𝑉_𝑒,𝑛+𝑉_𝑒,(𝑛+1))−𝐼𝑖,𝑛) (1)
where 𝑉_𝑒,𝑛 is the extracellular voltage at each node, 𝐶_𝑚 is the membrane capacity, 𝐺_𝑎 is the axoplasmic conductance, and 𝐼_𝑖,𝑛 is the ionic current across the membrane at each node. The ionic current was approximated as the sum of the current through 3 types of voltage gated ionic channels plus a leakage current as in (McIntyre et al 2002) (see Appendix A).
(...)
Equations (1) and (2) [equations for the end nodes] and the equations describing the evolution of the ionic currents were integrated by the implicit backward Euler method in MATLAB. The values of 𝑉_𝑒,𝑛 for each configuration were taken from the solution of the FEM model introduced in the previous section.

The Appendix contains exactly the same equations as the original paper by McIntyre et al., with (almost) exactly the same parameters.

So I thought that it should be possible to start with just a single node with those equations (contained in axnode.mod in the McIntyre model from the ModelDB #3810), which should remain at -80 mV if initialized properly and without a stimulus (which would correspond to zero voltage differences in equation (1)), and should fire with an apropriate stimulus.
Since this seems not to be the case, I must be missing something.

Thank you for helping me to find my error...

Matthias

[ted]

In the field of computational neuroscience, there is often a significant difference between the published narrative description of a model and the actual model that the authors used to generate their published results. That's one of the reasons why ModelDB exists.

You say your implementation of a model node of Ranvier is computationally equivalent to McIntyre et al.'s model node of Ranvier.
Mercadal et al. say their implementation is computationally equivalent to McIntyre et al.'s.
Which claim is correct?

You don't have Mercadal et al.'s source code so you can't test their claim. However, you can download the source code for the McIntyre et al. 2002 model and see what happens if stimulus amplitude is reduced to 0 nA. Is the McIntyre et al. 2002 model spontaneously active in the absence of stimulus?

[mzenker]

Hi Ted,

thank you for your message.
Perhaps I was not clear enough in my description of what I have done.

I have downloaded the McIntyre et al. source code from the ModelDB and run it with NEURON. If I don't change anything in the model, I see zero activity, i.e. a flat line at -80 mV when the stimulus is 0 nA. If I set the stimulus to 2 nA, it fires like it should.

Then I have built a node of Ranvier, using the same parameters (length, diameter) as McIntyre et al., and their axnode.mod which contains the modified HH equations, omitting the other parts of the model which describe the myelinated part of the fiber (see the hoc code in my first post). I did not change anything in the axnode.mod file. If I run that reduced model, I see permanent firing in absence of a stimulus.
(BTW before downloading NEURON, I have tried to solve those equations with Python. The result was essentially the same.)

As far as I can tell from their paper, Mercadel et al. claim to have done essentially the same, i.e. to use just the modified HH equations and not the myelinated part of the McIntyre et al. model, except that they use MATLAB instead of NEURON to solve the equations. I have sent them an e-mail asking for details some time ago, but did not receive an answer.

I see essentially two possible reasons why the node with the modified HH equations behaves the way it does in my implementation:

1. I have made mistakes or erroneous assumptions in the implementation.

2. If my implementation is correct, then those equations may be parametrized in a way that they don't work without the other parts of the model, meaning that the node model alone is not valid. That would mean that Mercadel et al. have done things they don't write in their paper to make it work for them.

My first guess would be reason no. 1. That's why I have posted my question here... ;)

If, however, reason no. 2 is the one to retain, then my next attempt would be to use the full model of McIntyre et al. and add extracellular stimulation, as described in this thread: viewtopic.php?f=8&t=1814.

My ultimate goal is to model the response of a motor neuron to an extracellular stimulus generated by two electrodes at some distance - just as in the paper by Mercadel et al., which is why I try to reproduce their model.

Thanks for comments...

Matthias

[ted]

As far as I can tell from their paper, Mercadel et al. claim to have done essentially the same, i.e. to use just the modified HH equations and not the myelinated part of the McIntyre et al. model

They represent internodes with resistors, which amounts to assuming that the effect of myelination is the same as surrounding the axon with a perfect insulator that has a vanishingly small dielectric constant.

I have sent them an e-mail asking for details some time ago, but did not receive an answer.

You might try again.

At this point, you don't know whether an isolated node of the McIntyre et al. 2002 model is spontaneously active or not. You have their source code, so why not do a simple computational experiment that will resolve this issue?
topology()
confirms that node[0] is the root section of the model (the section that has no parent), and that Its child is MYSA[0]. All you have to do is electrically decouple these sections from each other.

How?

node[0] psection()
reveals the properties of node[0], and
MYSA[0] psection()
reveals the properties of MYSA[0].
Both sections have nseg == 1, so you can make current flow through the cytoplasm negligible by setting their Ra parameters to 1e9 ohm cm. You'll also want to make extracellular current flow negligible, just for the sake of completeness.

Code: Select all

node[0] for i=0,1 print xraxial[i]

and

Code: Select all

MYSA[0] for i=0,1 print xraxial[i]

show that xraxial[1] is already 1e9 megohms/cm in both sections, so you only have to set their xraxial[0] parameters to the same value.

At this point you will have accomplished a virtual microdissection of a single node of Ranvier from the McIntyre et al. 2002 model axon, and you'll be ready to run the simulation that will reveal whether that isolated node is spontaneously active.

Suggestion for how to make the necessary parameter changes (assumes you have downloaded and expanded MRGaxon.zip, then compiled its mod file):
1. In the same directory, create a hoc file called changeparams.hoc that contains the statements that assign the desired parameter values for node[0] and MYSA[0].
2. Use NEURON to execute mosinit.hoc. This will create the model axon and set up the user interface.
3. At the oc> prompt execute
xopen("changeparams.hoc")
and now you're ready for the moment of truth.

[mzenker]

Hi Ted,

thank you very much for your helpful message. I have learned more about NEURON and about my case.
I have done what you suggested. I have created a file changeparams.hoc containing the following code:

Code: Select all

node[0] {
	xraxial[0]=1e9
	Ra = 1e9
	//Ra=rhoa/10000
	}
MYSA[0] {
	xraxial[0]=1e9
	Ra = 1e9
	//Ra=rhoa*(1/(paraD1/fiberD)^2)/10000
	}

Then I have started mosinit.hoc from the McIntyre et al. model and executed changeparams.hoc from the NEURON command line.
By playing with different combinations for Ra in the node and the MYSA sections, I find that Ra = 1e9 in either of them makes node[0] fire without stimulus, just as I have seen in my "node-only" model.

So the node model seems to need the connection to a neighboured section with finite cytoplasmatic resistivity in order to be in a steady state at the resting voltage. Thank you for helping me clarify this!
I am not sure if this was to be expected, however...

Now I could try to add an extracellular stimulus to the full McIntyre et al. model as described in the other tread, or build a simpler myelinated section using just a resistance, like Mercadel et al. seem to do it. I will do some further reading before I proceed.

BTW is there a trick to interrupt the output of topology(), or pipe it into a file? On my terminal (I am under Windows), the output is so long and large that I cannot scroll back to its beginning...

Thank you again for your help,
Matthias

[ted]

I am not sure if this was to be expected, however...

If we could predict everything, there would be no need to do any experiments of any kind (including experiments on computational models).

By playing with different combinations for Ra in the node and the MYSA sections, I find that Ra = 1e9 in either of them makes node[0] fire without stimulus, just as I have seen in my "node-only" model.

If you want to verify that your model is a correct implementation of McIntyre's model node of Ranvier, use the Vector class's record() method to capture the time course of each model's node.v(0.5), and compare the recorded values. The models differ in implementation, so results won't be absolutely identical, but it would be adequate to plot the two time courses and examine visually for significant differences over the first 10-20 ms.

So the node model seems to need the connection to a neighboured section with finite cytoplasmatic

cytoplasmic

resistivity in order to be in a steady state at the resting voltage.

The node model probably needs an electrical load in order to even have a resting potential. Notice that all of the non-node sections have not only pas but also membrane capacitance--both of which contribute to the stability of the intact axon model.

Take a close look at the spike train produced by the isolated node. Examine the first 15 ms closely, especially the interspike intervals. During the depolarizing phase of each ISI the net ionic membrane current i_ion must be inward (obviously). Not only is i_ion inward, but the more depolarized that v becomes, the faster it depolarizes, so the inward i_ion must be growing larger. After all, the node has voltage-dependent sodium channels that are activated by depolarization, and those account for acceleration of the depolarizing phase of the ISIs.

e_pas is -80 mV, so if you attach a section that has the pas mechanism to a node of Ranvier, that passive section will produce an outward current that tends to counteract the net inward current generated by the node. Furthermore, the more depolarized node.v becomes, the larger the passive section's i_pas will become, reducing the slope of the depolarizing phase of the ISIs and slowing the firing rate. If g_pas is big enough, i_pas will be too large for the node's inward current to drive the node to spike.

The passive section's membrane capacitance also reduces the rate of change of v (remember that dv/dt = (1/c)*net current), and this will have a flattening effect on the ISIs and slow the rate of spiking. If its membrane capacitance is big enough, dv/dt will be so slow that voltage-gated potassium channels have time to open and prevent the sodium channels from generating a spike.

You could explore the role of internodal cm and pas in the McIntyre model by a couple of simple computational experiments.
To make this easy, first create a set of sections that contains only internodal sections.

// make a set that contains all non-node sections
objref internodal
internodal = new SectionList()
forall internodal.append()
forsec "node" internodal.remove()

1. What happens if cm of all internodal sections is negligibly small?

// reduce cm of all non-nodal segments by a factor of 1e6
forsec internodal for (x,0) cm(x)*=1e-6

Now run a simulation; you might want to increase Tstop to 1500 or 2000 ms, just to be sure of stability.

2. What happens if g_pas of all internodal sections is negligibly small?

// restore cm to its original values
forsec internodal for (x,0) cm(x)*=1e6
// reduce g_pas by a factor of 1e6
forsec internodal for (x,0) g_pas(x)*=1e-6
// now run a simulation

You may want to increase Tstop to 1500 or 2000 ms.

3. What happens if both cm and g_pas of all internodal sections is negligibly small?

is there a trick to interrupt the output of topology(), or pipe it into a file? On my terminal (I am under Windows), the output is so long and large that I cannot scroll back to its beginning...

In recent versions of Windows, can't length of history be specified as one of the "preference" settings for a cmd window? What about NEURON's own bash terminal?

[mzenker]

Hi Ted,

thank you for yet another informative message and pedagogical challenge! ;)

I am not sure if this was to be expected, however...

If we could predict everything, there would be no need to do any experiments of any kind (including experiments on computational models).

Agreed. However, it is usually good to have an expectation for the outcome of an experiment or simulation. If the outcome doesn't match the expectation, then either the expectation is wrong (as it seems to be the case here), or there is an error in the experiment or simulation. One should never believe blindly in the correctness of results of experiments or simulations...
(Nature is always correct, of course, but experimental setups my contain errors.)

If you want to verify that your model is a correct implementation of McIntyre's model node of Ranvier, use the Vector class's record() method to capture the time course of each model's node.v(0.5), and compare the recorded values. The models differ in implementation, so results won't be absolutely identical, but it would be adequate to plot the two time courses and examine visually for significant differences over the first 10-20 ms.

That is what I had already done in the previous post. The result is that the firing patterns are visually identical between my reduced node model and node[0] of the McIntyre et al. model with Ra=1e9.
But I will finally use the original McIntyre model anyway.

So the node model seems to need the connection to a neighboured section with finite cytoplasmatic

cytoplasmic

Oops, sorry. I am german. And I am a physicist. ;)

resistivity in order to be in a steady state at the resting voltage.

The node model probably needs an electrical load in order to even have a resting potential. Notice that all of the non-node sections have not only pas but also membrane capacitance--both of which contribute to the stability of the intact axon model.
(...)
You could explore the role of internodal cm and pas in the McIntyre model by a couple of simple computational experiments.

I assume that those are meant be run on the original MRGaxon model with Ra set to its initial values, and without a stimulus. The observed quantity is node[10].v(0.5).

1. What happens if cm of all internodal sections is negligibly small?

No change. A flat line at -80 mV.

2. What happens if g_pas of all internodal sections is negligibly small?

v is increasing slowly and starts to exhibit spikes with a low repetition rate (delta t ~= 100 ms) at about 1200 ms.

3. What happens if both cm and g_pas of all internodal sections is negligibly small?

v exhibits quickly repeating spikes (delta t ~= 1 ms) from the beginning, in fact just as the isolated node.
So I have again learned something - thanks for explaining and helping me to gain more insight!
The bottom line for me is that the McIntyre et al. model is parametrized so that it does not make sense to use the node alone without also modeling the myelinated part in some way. Indeed, it is intended to represent a myelinated nerve containing nodes AND myelinated parts between the nodes. So my intention to get started by modeling the node alone was doomed to failure.

In recent versions of Windows, can't length of history be specified as one of the "preference" settings for a cmd window? What about NEURON's own bash terminal?

Hmmm, I didn't try the settings for the cmd window. But I could view the root of the topology by using the bash terminal, start the NEURON session from the command line with a "> logfile.txt" appended and execute the topology command blindly. Another way is to reduce the number of nodes to 11 (which is the minimum the model accepts). Then the topology plot fits better into the terminal. So problem solved. :)

I think the best bet for my purpose is finally to use the full McIntyre et al. model, which seems to be confirmed and reused by many researchers since more than 15 years, and just add an extracellular voltage stimulus. This has been discussed in several threads where you have kindly helped the posters, so I think I will succeed to get it to work. If not, I might come back with more questions.

Thank you again for your help and support, and for providing NEURON to the community!

Matthias

[ted]

Oops, sorry. I am german. And I am a physicist.

"Not that there's anything wrong with that," as Jerry Springer would say. OK, kytoplasmatisch.

A comment to those who might read this thread: How could the Matlab modelers escape problems with their own implementation, in which nodes were coupled by ordinary resistors (which omits representation of both membrane capacitance AND membrane leak current)? Whatever they did, they forgot to mention it in the article. This is yet another demonstration of why code sharing is important.