www.neuron.yale.edu

I have adapted a model of a myelinated axon (MRG axon - McIntyre et al, 2002) available on the modelDB so that it would include stimulation by both an extracellular point source and an intrinsic current clamp.

The current clamp will be injecting a bursting pattern on one end of the fiber and I want to record what happens on the other end of a population of fibers located at random distances from the electrode (extracellular point source) when stimulated at by a monophasic pulse at diferent frequencies.

The model behaves has expected for frequencies below 1kHz with the fibers closer to the electrode being driven by the extracellular stimulation and suppressing the bursting pattern.

For frequencies above that value (2kHz and 5kHz for instance) a strange phenomena happens. Fibers appear to be compleatly desynchronized and not following the extracellular stimulation at all. Some of those fibers could be going trough some sort of high frequency block similar to what can be seen in Kilgore et all 2006. However other appear to spike randomly. Please note that spikes appear always at the same time if I repeat the simulation with the same initial conditions.

I was expecting that with deterministic models of Na and K channels the fibers would fire with a frequency similar to 1/(node refractory period) for 2000Hz stimulation frequency.

By running the model with different timesteps I found that the higher the timestep the less of this pseudo-random desynchronization effect I could be seen.

On this page http://duke.edu/~jd139/neuronforum.html the rasterplots of how the same population of fibers behaves for high-frequencies and how that behavior varies with the time-step can be seen.

I am wondering if the effect that I am seeing is related to intrinsic model dynamics (cable equation and Hodgkin-Huxley equations) or whether it is related with the integration method and chosen timestep once I am using fast Na, slow K and persistent Na channels built on NMODL and solving the states with the "cnexp" method.

Thanks for your time, Joao.

On this page http://duke.edu/~jd139/neuronforum.html the rasterplots of how the same population of fibers behaves for high-frequencies and how that behavior varies with the time-step can be seen.

I am wondering if the effect that I am seeing is related to intrinsic model dynamics (cable equation and Hodgkin-Huxley equations) or whether it is related with the integration method and chosen timestep

The following

By running the model with different timesteps I found that the higher the timestep the less of this pseudo-random desynchronization effect I could be seen.

certainly indicates that the answer to your question is "numerical methods account for at least part of what you saw." dt 0.025 ms is barely short enough for good accuracy in simulations that involve the HH mechanism, which has slower dynamics than axonal ion channels in vertebrates. If you see unaccountable firing patters, try an even smaller dt. Does your model implementation work with cvode on, or with secondorder=2? Either of those options might give you increased numerical precision without having to reduce dt to excessively tiny values.

Two items that I forgot to mention:

First, the title of this thread is incorrect. The spike times might be chaotic, or appear to be chaotic, but they are not "random" because they are generated by a completely deterministic system--as evidenced by the fact that, given a particular set of initial conditions, results are reproducible from run to run.

Second, it is known that (completely deterministic) nonlinear systems driven by periodic inputs can produce chaotic output--a topic of considerable practical importance, especially with regard to cardiac arrhythmia and pacemakers.