www.neuron.yale.edu

Hi,

I've been simulating various transcranial magnetic stimulation (TMS) waveforms for extracellular stimulation in NEURON using the xtra.mod file approach written by Ted. Typically, I down-sample a recorded waveform to the time step of my simulation (e.g. 1 µs or 5 µs) in Matlab, and then I use vector.play to "play" the normalized waveform into the global "stim" variable with continuous mode on:
I noticed that if, for example, I use a TMS waveform recorded at 1 µs steps and a NEURON simulation time step dt of 5 µs for either of the two following cases, I get different results:
1) Downsample TMS waveform in MATLAB to 5 µs and import to NEURON (downsample.m function, just takes every 5th sample in this case), or
2) Import original TMS waveform (1 µs), and let NEURON handle downsampling at simulation time steps

And by different results, I mean different polarizations for a given compartment in the neuron model at the same stimulus amplitude. Is this due to NEURON using the additional waveform samples between time steps? And if so, is it better practice to just import the recorded waveform with the highest temporal resolution available, rather than downsampling to the simulation time step?

The other concern I had was avoiding anti-aliasing of the waveform, so I apply a 1st order, 200 kHz low-pass filter in MATLAB to reduce anti-aliasing from the subsequent downsampling for simulations. Does NEURON do anything under the hood to deal with aliasing?

if, for example, I use a TMS waveform recorded at 1 µs steps and a NEURON simulation time step dt of 5 µs for either of the two following cases, I get different results

How different is "different"? Major qualitative difference, or just a noticeable quantitative difference while overall time course is preserved?

Is this due to NEURON using the additional waveform samples between time steps?

I initially misread this question. In a simulation the driven variable has a constant value during each fadvance(). If dt, the integration time step, is the same as the data sampling interval dts (and no roundoff error issues arise), then the driven variable's values will be x0, x1, x2, x3, x4, x5, . . . If dt is longer than dts, some sampled values will be skipped. For example, if dt is exactly equal to 5*dts, the driven variable's values will be x0, x5, x10 etc.. At the very least, you can expect that the high frequency components of the model's response will be altered. The lower frequency components are also likely to be different; after all, it is unlikely that x5 will even equal the mean value of x1..x4.

There might also be some "finite precision" issues if the sampling times cannot be expressed exactly with double precision arithmetic, so that the driven variable's values are occasionally at 4 or 6 sample intervals.

Do the recorded data consist of a single stream of values (samples of the measured variable), or two streams (one being samples of the measured variable and the other being samples of the corresponding time)? At what point in the simulation do the differences first appear? Within the first 1000 samples, or after many thousands or tens of thousands of samples? If the recorded data include sampling times are you using those times to drive Vector.play? Do you get better results if you use NEURON to compute the tvec that is supplied as an argument to Vector.play?

is it better practice to just import the recorded waveform with the highest temporal resolution available, rather than downsampling to the simulation time step?

Using the "highest temporal resolution available" allows roundoff error to determine which samples are fed to the model. It's probably better to start with data captured at intervals that are an integer fraction 1/N of dt, and feed every Nth value to NEURON. Or, if necessary, used a principled algorithm (even linear interpolation if necessary) to resample at intervals of dt.

The other concern I had was avoiding anti-aliasing of the waveform, so I apply a 1st order, 200 kHz low-pass filter in MATLAB to reduce anti-aliasing from the subsequent downsampling for simulations. Does NEURON do anything under the hood to deal with aliasing?

Nope, but there's probably more than one Python library that has methods for antialiasing time series data.

How different is "different"? Major qualitative difference, or just a noticeable quantitative difference while overall time course is preserved?

More the latter, probably <5% difference in most cases I looked at.

Do the recorded data consist of a single stream of values (samples of the measured variable), or two streams (one being samples of the measured variable and the other being samples of the corresponding time)?

2nd one. I input a time vector and Efield vector. In the 1st case (downsample in MATLAB), this should be the same time points as the simulation, and in the 2nd case (import original waveform), the time vector is at 1 µs steps, while the neuron is simulated at 5 µs steps. The time values should should be identical to a time vector generated in NEURON, since it's at regularly spaced intervals (generated using colon operator in MATLAB), not using output from the oscilloscope.

At what point in the simulation do the differences first appear?

Within the first few samples. I plotted voltage recordings for the cases I mentioned (1st image below), as well as for simulations with dt = 4 µs, so I could easily downsample to 4 µs and 2 µs (2nd image below). It looks like inputting the waveform downsampled to 2 µs produces the same membrane voltages as inputting the original waveform at 1 µs, so I think that explains what's going on. I'm thinking I can just make sure to sample my waveform with steps <dt/2 and input that to my NEURON simulation, or like in this case, just keep it at the 1 µs temporal resolution, since it seems to give the same result.

NEURON's defaut integrator is implicit Euler, which has first order accuracy. This means that a driven variable, be it a stimulus amplitude or a time-varying parameter, is effectively a staircase waveform--changes abruptly at each time step boundary, and remains constant until the next time step boundary. Consequently extremely fast transitions will be closely approximated only if dt is 5-10 times shorter, and fast oscillations will be poorly preserved unless dt is shorter than ~0.1 * the period.

Fortunately, neurons have time constants on the order of 10s of ms, so all that superfast stuff is profoundly smoothed, as your simulated membrane potential shows.

I'm thinking I can just make sure to sample my waveform with steps <dt/2 and input that to my NEURON simulation, or like in this case, just keep it at the 1 µs temporal resolution, since it seems to give the same result.

That's a sensible conclusion. Addressing these kinds of questions often requires such an empirical approach. After all, we're trying to approximate a system that is continuous in space and time with a bagful of equations that sample time and space at a finite number of points. Results will always be approximate, and interpretation and judgement must always be applied to decide what is "good enough". The acceptable goal is robustness of important qualitative outcomes in the face of perturbation of model and simulation parameters.