A general method to generate artificial spike train populations matching recorded neurons

Abstract

We developed a general method to generate populations of artificial spike trains (ASTs) that match the statistics of recorded neurons. The method is based on computing a Gaussian local rate function of the recorded spike trains, which results in rate templates from which ASTs are drawn as gamma distributed processes with a refractory period. Multiple instances of spike trains can be sampled from the same rate templates. Importantly, we can manipulate rate-covariances between spike trains by performing simple algorithmic transformations on the rate templates, such as filtering or amplifying specific frequency bands, and adding behavior related rate modulations. The method was examined for accuracy and limitations using surrogate data such as sine wave rate templates, and was then verified for recorded spike trains from cerebellum and cerebral cortex. We found that ASTs generated with this method can closely follow the firing rate and local as well as global spike time variance and power spectrum. The method is primarily intended to generate well-controlled spike train populations as inputs for dynamic clamp studies or biophysically realistic multicompartmental models. Such inputs are essential to study detailed properties of synaptic integration with well-controlled input patterns that mimic the in vivo situation while allowing manipulation of input rate covariances at different time scales.

Keywords: Cerebellum; Correlation; Cortex; In vivo; Model; Synaptic.

Fig. 1

AST targets. a) Adaptive Gaussian Rate template from a 10s Purkinje cell spike train (see Fig. 2 for construction process). b) Interspike-Interval (ISI) histogram of Purkinje cell recording (7474 spikes see Table 1). Spike density is calculated as spikes per bin per second. c) Power spectrum of Purkinje cell spike train.

Fig. 2.

Construction of adaptive Gaussian Rate Template a) 5s segment of spike train from area MT neuron, b) Each spike (dotted vertical lines) is first convolved with a fixed width Gaussian and the Gaussians are summed up (green line, slow rate function). Then the slow rate template is used to determine an adapted width of a Gaussian for each spike to produce the adaptive rate template (see Methods for details).

Fig. 3.

Examples of rate template manipulations to control AST properties. a) aGLR of a 5 s segment from a Purkinje cell recording (black) and a manipulated rate template (red) in which the 3–5 Hz frequency component has been reduced by applying a gain of 0.1 to this frequency band (for details see our Matlab function splitsignal_adjust() at https://doi.org/10.15139/S3/8ILYHZ. b) Power spectrum of the original rate template and the 2–5 Hz filtered template. Data are from a 115 s template. c) aGLR from a 10 s segment of the same Purkinje cell recording (black) and a manipulated rate template (red), in which a beta frequency band (12 −16 Hz) was amplified by a gain factor of 20. d) Power spectrum of beta-band amplified rate template. Data are from a 115 s template.

Fig. 4.

Matching recorded spike trains to ASTs for PC, MF and MT neurons. a-c) For each cell type a segment of spike data is shown from the recording (REC, top trace) and an AST generated from the aGLR of the recorded spike train. Note that the time base is different for each neuron based on their mean spike rate. d-f) The aGLR is shown for a segment of the full template for the recording of each cell type (black trace), and superposed the average aGLR from 100 ASTs is depicted (red trace). g-h) The ISI distribution of the full length recordings (PC: 115s, MF: 45s, MT42: 1,079s) and the mean ISI distribution of 100 ASTs of the same length (red). For these simulations, physiologically plausible assumptions were made about the free parameters of our search algorithm. See table 4 for error comparison to other parameter settings.

Fig. 5.

Manipulating rate co-correlations within populations of ASTs. a) The black trace shows the ‘slow’ rate template with a σ of 100 ms from the 115s PC recording. The red traces show the average rate template also computed with a 100 ms σ from 100 ASTs where each AST was drawn from the same PC aGLR. b) Same black trace, but 50 of the 100 ASTs for which the mean ‘slow’ rate template is shown (red) are drawn from a time-shifted version of the PC aGLR. c) All 100 ASTs are drawn from time shifted aGLRs. d-f) Average sliding window peak correlation at lag 0 between ‘slow’ rate templates computed for each AST and the rate template from the recording. Note that the correlations shown here are lower than they appear in panels a-c above because individual ASTs carry substantial individual fluctuations and are less correlated to the recording than the average rate template of all ASTs where these individual random fluctuations cancel out. Note that the sliding window correlation shows lower values at times when the rate change amplitude in the recording are lower as the noise-driven fluctuations in the AST rate at those times are relatively more dominant. The average peak correlation for an SF=0 was 0.70, for SF=0.5 was 0.41, and for SF=1 was 0.05.

Fig. 6.

Parameter search for best algorithm performance. abc1–3) For each cell type we plotted heat maps of the average percent error of 100 ASTs against the original recorded spike train for firing rate (FR), CV, and the average per bin deviation in the ISI density over 110 (PC), 150 (MF) or 200 (MT42) 1ms bins starting at t=0 ms ISIs (see Fig. 4). Note that the error ranges are quite different for different plots, and the color bar ranges are separately adjusted from min to max error for each plot. abc4) The standard deviation of the total error (rate + CV + LV + ISI errors) is shown for different population sizes of ASTs for which these properties were averaged. AST populations of size were constructed 100 times, and the standard deviation between populations was calculated. As the total error was calculated as a % deviation from the mean of the recorded measures, a value for example of ‘2’in the StDev measurement denotes that the total error had a standard deviation of 2%. See table 3 for the absolute size of the different errors.

Fig. 7.

Constant target rate templates. Surrogate data were constructed as rate templates with a constant target rate and low or high spiking regularity (LV of 0.1 or 1.5). We find that the average spike rate from 100 ASTs matches the target rate with only very small errors (see Table 2). The individual ASTs (blue traces) show large rate fluctuations due to the stochastic nature of the gamma spike train. These random rate fluctuations are much more pronounced for a highly irregular spike train target (Fig. 7b,d; LV=1.5) vs. a regular spike train target (Fig. 7a,c). Individual AST rates were constructed as aGLRs with a scale factor (sf, eq. 4) of 1.0 instead of 0.25 as used for physiological spike trains in order to more clearly show the contribution of single spikes to rate changes (Fig 4a,b; blue trace).

Fig. 8.

Sinusoidal target rate templates with increasing frequency (zap). Black traces: Surrogate data were constructed as rate templates with a sinusoidal rate change and increasing frequency, and given a low or high target regularity (LV of 0.1 or 1.5). a-d) Red traces: Average aGLR of 100 ASTs computed from the target template at low sinusoidal frequencies. A sample AST is depicted with black bars at the bottom of each panel. e-f) Sinusoidal template (black) and average AST aGLR for higher sinusoidal frequencies and target spike rates of 10–20 Hz (trough and peak of sinusoidal). g) Same for a target spike rate of 40–100 Hz (trough and peak of sinusoidal) and a regular spike train target (LV = 0.1).

Fig. 9.

Physiological spike train templates with superposed rate change events. Rate templates from our 3 sample neuron types were convolved with a square rate change event of 0.4s (a-c) or 2s (d-f) duration. Square rate changes were +− 0.5* mean spike rate for a modulation strength of 1 (a-c, e). The rate change was applied to the template at even intervals of 5s for the 2s duration rate changes and every 1s for the 400ms rate changes. 100 ASTs were constructed from the recorded spike train with the superposed square waveform rate changes, which represent a canonical form of a behavioral event related spike rate change. We then constructed peri-event time histograms (PSTH) aligned to the rate changes. a-f). The top panel shows a raster plot of 50 spike trains aligned to the center of the square rate change. The dot size shown for each spike is scaled to the mean spike rate. The bottom panels show the average rate template around the imposed spike rate change (blue traces), which is the target rate for the ASTs. Black trace: The mean spike rate per 1ms bin of all event aligned spike trains (N=11,200 [2100] for PC, N=4,400 [900] for MF, and N=107,900 [21,600] for MT42), where [] brackets denote the event numbers for panels d-f. Red traces: Mean event aligned aGLR constructed from the ASTs.