Cross Nearest-spike Interval Based Method to Measure Synchrony Dynamics

A new synchrony index for neural activity is defined in this paper. The method is able to measure synchrony dynamics in low firing rate scenarios. It is based on the computation of the time intervals between nearest spikes of two given spike trains. Generalized additive models are proposed for the syn-chrony profiles obtained by this method. Two hypothesis tests are proposed to assess for differences in the level of synchronization in a real data example. Bootstrap methods are used to calibrate the distribution of the tests. Also, the expected synchrony due to chance is computed analytically and by simulation to assess for actual synchronization. 1. Introduction. Neurons carry information through the nervous system by means of Action Potentials (AP), also called spikes. Action Potentials are electrical pulses of similar amplitude which last approximately for 1 ms. Spike trains are sequences of such AP. Current models of brain processing indicate that neuronal information is mostly coded in the time course of AP, i.e., the firing rates, and in the temporal associations between spike trains of different neurons. Temporal correlation and synchrony are therefore key features to understand neuronal coding. Depending on the state of the brain, the firing rate may change considerably. Usual studies deal with neuronal data with high activity, as the individuals are performing tasks. But, high firing rates are not always the case. For example, during slow-sleep, neurons present a low firing rate, usually grouped in bursts with a frequency of 1-5 Hz (the so called delta rhythm). On the contrary, in the awake state, or while the individual is performing a task, firing rates increase and the AP responses turn to a characteristic tonic mode. Synchrony in neuronal activity corresponds to the joint activity of neurons, that is, neurons firing close together in time. Electrophysiological studies commonly use

1. Introduction.Neurons carry information through the nervous system by means of Action Potentials (AP), also called spikes.Action Potentials are electrical pulses of similar amplitude which last approximately for 1 ms.Spike trains are sequences of such AP.Current models of brain processing indicate that neuronal information is mostly coded in the time course of AP, i.e., the firing rates, and in the temporal associations between spike trains of different neurons.Temporal correlation and synchrony are therefore key features to understand neuronal coding.
Depending on the state of the brain, the firing rate may change considerably.Usual studies deal with neuronal data with high activity, as the individuals are performing tasks.But, high firing rates are not always the case.For example, during slow-sleep, neurons present a low firing rate, usually grouped in bursts with a frequency of 1-5 Hz (the so called delta rhythm).On the contrary, in the awake state, or while the individual is performing a task, firing rates increase and the AP responses turn to a characteristic tonic mode.
Synchrony in neuronal activity corresponds to the joint activity of neurons, that is, neurons firing close together in time.Electrophysiological studies commonly use 28 GONZ ÁLEZ, CAO, FAES, MOLENBERGHS, ESPINOSA, CUDEIRO AND MARI ÑO synchrony measures to reveal computing properties among neurons under different conditions, such as different kinds of stimuli.In this work we present a synchrony measure that can be easily adapted to different firing rate scenarios.The objective of the experimental work from which we take the data is to study synchrony dynamics under low firing rates and to distinguish differences in synchrony profiles.Hence, we show the performance of our method measuring synchrony in real data, recorded in a repose state of spontaneous oscillatory activity which is briefly disrupted by the experimental activation of two neural pathways.We discuss tests to determine (1) whether there is a significant change in synchrony due to electrophysiological stimulation and (2) whether the synchrony observed under the effect of the stimulation differs between experimental conditions.Also, we compare the observed amount of synchrony with the amount of synchrony expected due to chance when there is firing activity, to assess whether the observed synchrony is relevant.On the other hand, there are experimental preparations, such as some electrophysiological recordings "in vivo", in which getting many recordings of one group of neurons may be very difficult and hard to attain.The proposed method can be used to measure synchrony in both high and low firing rates, and with small numbers of trials (experimental repetitions).We also examine the performance of the method by means of a simulation study.
There have been proposed numerous methods devoted to measure synchrony, many of them are based on cross-correlation analysis.Most cross-correlation indexes can be applied to neuronal spontaneous activity (typically with a low firing rate), provided that the recording period is long enough to get a minimum amount of spikes, but those tools are not useful to calculate instantaneous synchronization dynamics, because of the characteristic low amount of spontaneous spike activity.For example, the joint peristimulus time histogram (JPSTH) [6] displays the dynamics of correlation between neurons.Its normalized version is the Pearson correlation coefficient (computed across trials) of the firing counts of both neurons at two different time bins.Other methods commonly used to capture synchrony are those based on 'unitary events' [7,8].These methods rely on binned trains, where, for each, bin activity is denoted by a 1 if a spike occurs in the interval or 0 otherwise.Unitary events refer to the occurrence of coincident spikes, or 1-1 matches, in the neurons under study.The unitary events analysis estimates the probabilities of joint-firing under the hypothesis of independence of the two spike trains.These probabilities are used to compute the expected numbers of joint spikes.The test for presence of synchrony is defined in terms of the difference between the expected frequencies and the observed ones.Faes and co-authors, in [5], propose a synchrony index, the Conditional Synchrony Measure, which is calculated, also, with binned trains.It is a flexible method, based on estimating the probability of joint-firing given that one of the neurons fired.However, it is not useful in low firing rate scenarios.Another method is presented by Quiroga et al. in [13].This method is based on relative timings of events in time series.For general spike train analysis methods and state-of-the-art, including correlation among neurons, see [3] and [10].
2. Data.We illustrate the performance of our method with a real-data example.These data corresponds to a low firing rate scenario.In the experiments, simultaneous neuronal spike activity in the primary visual cortex of anesthetized cats was recorded.The anesthetized state attained in our experiments is very similar to the slow-sleep state and therefore neurons present an oscillatory low firing rate.The neuronal pattern of activity during sleep mode was artificially disrupted by the use of micro-electrical stimulation in two different regions of the brain: brainstem (bs) and basal forebrain (bf ).These two stimulated areas are known to modulate the sleep-wake cycle and are essential in the 'waking up' process affecting the whole brain, inducing a transition to a neuronal tonic state [1,11,14,15]).Three trials were recorded under each condition.The aim of the experimental study is to investigate the impact of bs and bf stimulation on the synchronization of neurons in the primary visual cortex.We will use the simultaneous activity recorded from two neurons under the just described experimental setting to illustrate our method.The length of the recordings considered in this study is 150 seconds, the stimulus onset occurring at 50 s.Although, the actual stimulation lasts for 2 seconds, its activating effects last longer, for up to 20 s. Figure 1 shows the firing rates of two trials (one under each condition) of the two neurons that will be used as an example.The firing rates, λ(t) were estimated using nonparametric kernel methods.The kernel estimator, λ(t), has the form: , where K is a kernel function and h is the bandwidth that controls the amount of smoothing [16].For illustrative purposes, we used a Gaussian kernel function with an ad hoc choice of the bandwidth.A study on the effective choice of the bandwidth can be found in [12].To correct for boundary bias we have used a mirror method, i.e., the data close to the boundaries were repeated in reverse order after the endpoint, prior to applying the smoothing procedure.In the top panel we can observe the firing rates of the two chosen neurons during a trial in which bs was stimulated (at 50 s).The bottom panel shows another trial in which bf was stimulated.In these four cases we can observe rates not higher than 3 Hz.
Figure 2 shows the JPSTH computed across the three trials from the bs stimulation, using two different choices of bin widths (1 ms: top panel; 100 ms: middle panel) and the synchrony measure proposed in this paper (bottom panel), which will be presented in Section 3. In the JPSTH analysis, the most commonly used bin width is 1 ms.If this is the case, there can only be one spike per bin because of the refractory period of neurons.Nevertheless, wider bins can be used and multiple spikes could be reduced to one.In the present example, when 1 ms bins are used, no information can be drawn from the JPSTH.When wider bins are used, some features start to be seen, such as larger coincidence rate right after the stimulus.Due to the small number of spikes and trials, information from wider time periods should be considered.However, when the JPSTH is the tool, increasing the bin width results in information loss, because the only important fact to compute the JPSTH is the presence/absence of spikes.In this context, when low firing rates and small number of trial induce to the use of larger time bins, the available information of the number of spikes that fall in each time bin is almost certainly valuable.This is why a method that takes into account the extra information is needed.Moreover, we consider that results will be more conclusive with such a method and hypothesis tests will be more capable of detecting differences in the association profiles.
3. CNSI based synchrony measure, CSM.In neurophysiological studies, two (or more) neurons are considered to be synchronized if there exists a close relationship between the time-course of their spike activity.Because of the broad range of neuronal activation, it does not exist a unique measure of synchronization, and experimental neurophysiologist adjust the acceptable conditions depending on the Firing rate (Hz) Firing rate (Hz) experimental conditions and the hypothesis under study.For example, the synchronization required to demonstrate a direct connection between two cells may require a time-lag between spikes of no more than 1 ms; but, on the other hand, during sleep oscillations, neurons are considered to be highly synchronized even with time-lags of tens of ms.During spontaneous activity, the overall activity in the brain, as can be observed in an EEG, presents very highly synchronous activity with large amplitude and low frequency waves, but, on the other hand, firing rates of individual neurons in this context are very low.Some statistical tools to measure synchrony fail in such a context because they are based on 1-1 matches, [5,7,8].This is, data is discretized into (usually) 1 ms bins, where the value at the i-th bin will be 1 if there is a spike and 0 otherwise.So, there is a 1-1 match between two spike trains at the i-th bin if there are spikes at that bin from both spike trains.In the context of low firing rates, the definition of synchrony should be relaxed by using wider windows than 1 ms.Let X = {X i } J1 i=1 and Y = {Y j } J2 j=1 be the spike times of two simultaneously recorded spike trains in the time interval [0, T ) and let {N X (t), t ∈ [0, T )} and {N Y (t), t ∈ [0, T )} be the counting processes associated to them.This means, N X (t) = #{X i < t, i = 1, 2, . . .J 1 } and similarly for N Y (t).Also, let n = J 1 + J 2 be the total amount of spikes in the two neurons.
Let us define the Cross Nearest-Spike Interval (CNSI) as the time that elapses between every spike of one neuron and the closest spike of the other neuron, not necessarily forward.Formally, let U −1 and U 1 be the waiting times from one spike As already discussed, when the firing rate of spike trains is very low, we need a broader definition of synchrony because exact firing matches hardly ever occur in 1 ms windows, not even in 10 ms windows.Given a pair of spike trains, we define synchrony as the event of two AP (one of each neuron) occurring within a time window of width δ.Note that δ is an intrinsic parameter to the measure and, therefore, needs to be chosen according to the biological problem under study.Nevertheless, issues about the selection of δ are discussed in simulated examples in Section 9. Let n δ be the number of CNSIs smaller than or equal to δ and define the following measure, which will be called CNSI-Synchrony Measure (CSM) and will be denoted p δ : This is a global measure of synchrony.It is symmetric and does not pick up causality.In terms of the counting processes, n δ can be expressed as: where I(A) is the indicator function of the event A.
In order to define CSM as a function of time, p δ (t), to take into account the non stationarity of the processes we define n δ (t) as follows: In practice, X j = t and Y j = t are events of probability zero but, even when working with binned trains, observing a spike in a given bin is very difficult when firing rates are low.So, in order to solve this problem, we will work with the information provided by a neighborhood of t.Let V t = (t − v, t + v] be a symmetric time window of length 2v around t, with v chosen by the researcher depending on the context.For higher firing rates, smaller windows could be used.Therefore, and For each t, we compute the CSM in the time window V t to obtain p δ (t).Consequently, to measure synchrony at a time point t 0 , the probabilities P (N Y (X j + δ) − N Y (X j − δ) ≥ 1) and P (N X (Y j + δ) − N X (Y j − δ) ≥ 1) are considered constant within V t0 and estimated from information in the whole interval.In this way, the local low firing rate is outweighed by neighborhood activity.
Population-wise talking, p δ (t) is an estimator of the probability, π δ (t), of (given two trains) finding two spikes (one of each neuron) closer than the quantity δ, given that occurred one spike at time t.In general, π δ can be considered as the probability of success of a certain Bernoulli trial, where the trial corresponds with one observation of the U with success probability that U is smaller than δ.Therefore, n δ becomes the observation of a Binomial variable, η δ , that counts the number of successes of the Bernoulli trials, η δ ∼ Binomial(N, π δ ) and consequently, p δ is an estimator of the probability π δ , p δ = πδ .Finally, the same argument can be used at each time window V t to obtain p δ (t) as the estimator of π δ (t) where η δ (t) ∼ Binomial(N (t), π δ (t)).In the rest of the paper we will drop the subscript δ except for the quantity n δ , because there is no possibility of confusion as everything depends on this quantity in the same way.
4. Model formulation.In this Section we present a model for the CSM.As already discussed, at each time window, the proposed measure is an estimator of the probability of a Binomial process where success is the event of observing Ũ smaller than δ and the total number of CNSIs is the number of Bernoulli trials, η(t) ∼ Binomial(n(t), π(t)).We propose to use Binomial Generalized Additive Models (GAM) with a logit link to explain the proportion of 'small' CNSIs in time.The convenience of using GAMs lies in the flexibility of non-parametric functions as it would be not reasonable to choose parametric response function.Also, the additive terms of the model allows for an easy interpretation.
In general GAMs can be represented as follows: and Y is a random variable belonging to the Exponential family of distributions, g is called a link function, S * i is a row of the model matrix for the parametric component of the model, α is the corresponding parameter vector and the f j are smooth functions of the covariates x j [9,17].
For our problem, we propose the use of binomial GAMs.In this particular case the canonical link function is used, that is, the logit: logit(µ) = µ 1−µ .Also, it is important to recall that we have several experimental conditions (stimuli) which are assumed to have different probabilities of joint firing.The model we propose is a mean curve plus difference one: where π j (t) stands for the success probability in η j (t) ∼ Binomial(N j (t), π j (t)) under the j-th experimental condition.The parametric part is reduced only to an intercept and f 0 (t) is a common curve for all the conditions, while f j (t) represents the difference from the curve f 0 (t) to the one corresponding to condition j = 1, . . ., J − 1.Finally, the t are the errors of the model which will be discussed later on.
In practice, we only have two conditions, J = 2, and we have three repetitions for each condition.Let n δji (t) be the observed n δ in the time window V t under the i-th trial of the j-th condition and n ji (t) the total amount of CNSIs in V t , j = 1, 2 and i = 1, 2, 3, and therefore p ji (t) = where, after obtaining the estimates, β0 , β1 , f0 and f1 we get the estimate πj (t) for each condition.We are interested in the time evolution of synchrony and therefore in the differences that could occur in certain periods of time.Nonetheless, a first overall study of the synchrony profile is cautious.For this matter, we will fit a simple model, a single curve model (Model 2), and compare it with Model 1: Model 2 in (3) corresponds to the hypothesis of no differences between experimental conditions over the whole time period.It comprises a single function for both conditions, while Model 1 in (2) allows the smooth functions to differ.With the present approach, in a two experimental conditions model, a possible difference between the two conditions is easier to observe because, if there are no differences between conditions then f1 (t) will be essentially flat.
To represent the smooth functions we have used penalized regression splines, with cubic splines basis.Smoothing parameters are chosen by the Generalized Cross Validation (GCV) criterion and the number of knots using the Akaike Information Criterion (AIC) value.The mgcv R package was used for this aim.
An important issue to note is that no dependency between time points is taken into account with these proposed models.Nevertheless, we are aware of the high dependency that exists between consecutive time points.The dependence will be taken into consideration when building confidence bands using the estimators and when developing the hypothesis tests.For this aim, we consider that the dependence is gathered in the errors and we propose an autoregressive model of order one, AR(1) for them: We can estimate γ by ordinary least squares.The estimated errors are ˆ ji (t) = p ji (t) − πj (t) and let γ be the minimizer of the sum of squares: which results in: ji (t − 1) and also we can estimate σ 2 by σ2 , the variance of the values âji (t) = ˆ ji (t) − γˆ ji (t − 1).

Hypothesis testing.
There are relevant questions about the synchrony between neurons to be asked.In most experiments, the different conditions refer to the application of a stimulus, the beginning of a task, or other kinds of event that have a starting point in time, which does not necessarily correspond to the beginning of the recording.If this is the case, it is of interest to investigate if synchrony is altered by the apparition of the aforementioned event.A second question to be asked is the possible differential effect of each condition on synchrony.Therefore, considering two experimental conditions, 1 and 2, the two main null hypotheses to test are -H 1 0 : π 1 (t) = π pre for every t after the condition onset.Here, π pre represents the baseline synchrony before the stimulus.A similar hypothesis can be formulated for condition 2.
-H 2 0 : π 1 (t) = π 2 (t).With this test we aim to detect differences in the synchrony profile induced by different experimental conditions.To test these hypotheses we will use the CSM as the test statistic and bootstrap critical bands to find out whether there are significant differences.The bootstrap procedure will be discussed in Section 7.
6. Synchrony due to firing rate.Another important aspect to investigate is whether the observed synchrony is different from the one that would be expected just by chance.Even in the case of two independent spike trains, CSM would not be zero, and its value could be large due to random close firing.It seems reasonable that this synchrony obtained by chance will increase with increasing firing rates.In the way the CSM is built, it could be affected by the firing rate and by its changes in time.To study this aspect we will compare the observed CSM with the one expected only by chance, due to the firing rate.For this purpose we need to calculate what CSM values are expected for the observed activity if the two trains were independent.These values will be approximated in this section in two different ways: 1) approximating numerically the theoretical expected synchrony by chance and 2) by simulating independent trains and modeling their synchrony.

Let us start with the theoretical expression for the CSM. Let
as an estimator of the expected value of: Taking expectations and using the property E(A) = E(E(A|B)) we have: which, under the assumption of independence, becomes: where ρ δY (s) = E(I{A Y (t, δ) = ∅}) which is the expected value for the indicator that there is a 'jump' in process N Y in (s−δ, s+δ].As we are using the information of the whole interval (t − v, t + v] to compute n δ , we can assume the process is stationary in this interval, so, ρ δY (s) ≡ ρ δY constant in (t − v, t + v].If this is not the case, but v is sufficiently small, then, ρ δY (s) can be approximated by a constant value, for example, the value of ρ δY (s) at the middle point of the interval: ρ δY (s) ≈ ρ δY (t).
From the previous paragraph we have, The previous discussion also holds for the quantity in (1), so finally we get that, under independence, where, On the other hand, we have that the total number of CNSIs, n, in the time window is an approximation to the CSM under the independence hypothesis.Note that the expected random synchrony under independence, given the firing rates of the two spike trains, can also be approximated using simulation.For this, we assume there exists an underlying function, g, of synchrony, dependent on the firing rates of the two spike trains, r 1 and r 2 , i.e., g(r 1 , r 2 ) = E (synchrony between trains Xand Y under independence|r 1 , r 2 ) .
Given the instant probability of firing, p 1 and p 2 , we simulate independent spike trains with the corresponding firing rates and compute the CSM for the simulated trains.This procedure is repeated for every pair (p 1 , p 2 ) in a two dimensional grid, P , spanning from 1 Hz to 100 Hz.We take steps of size 1 Hz for firing rates between every 2 Hz from 12 to 60 and every 5 Hz from 65 to 100 Hz.This function is then smoothed with a two-dimensional spline, and denoted by g(p 1 , p 2 ).
To simulate a spike train with a global firing rate of p 1 , 500 s time intervals are divided into 1 ms bins.In each bin, a Bernoulli trial with success probability p 1 (spiking) is drawn.Let {b i } M i=1 be the indices of the bins in which spikes are randomly assigned.Then, the spike times (in ms) for the simulated train are X * j = 0.001b j j = 1, . . ., M .The same procedure is used to simulate trains with global firing rate equal to p 2 .Using the two simulated spike trains with firing probabilities p 1 and p 2 , the CSM is calculated and an average over the 500 s is considered as the approximation of g at (p 1 , p 2 ).Notice that this simulation scheme is not independent of δ and v and therefore it needs to be carried out for each choice of these two parameters.
Finally, the smoothed trend is used to calculate the expected synchrony, under the assumption of independence, for the real pair of neurons.At each time point, the instant firing rate of each real train is estimated and those values evaluated in g(r 1 , r 2 ) to obtain the expected value.
7. Bootstrap confidence bands and testing for differences.A bootstrap procedure is carried out to build confidence bands for the predictions of the selected model.It is also used to build critical regions for the predictors under the null hypotheses described in Section 5, to assess for differences in synchrony.The bootstrap procedure is now described to build confidence bands for the estimators and later the procedure is slightly changed for hypotheses testing.The procedure for I trials is the following: 1. Fit the penalized regression model to the response data, p ji (t), and obtain the fitted probabilities, πj (t), for each condition j.

Fit the regression model from
Step 1 to the bootstrap data to obtain the bootstrap synchrony curve π * j (t). 7. Repeat Steps 4-6 M times to obtain M bootstrap curves for each condition: π * 1 j (t), . . ., π * M j (t).To build the confidence bands for the estimators, the only step left is to compute the α/2 and 1 − α/2 quantiles at each time point t to obtain (1 − α)% confidence bands at each condition j = 1, 2.
For the hypothesis tests, the procedure changes slightly mainly because the model to be fitted has to be the one under the null hypothesis.It is worth to note that the errors used to fit the AR model that will be used to build the bootstrap samples are the ones obtained by fitting the general model (2).Therefore, the procedure for I trials is the following: To fit the AR(1) for the errors: 1. Fit the general penalized regression model (as before) to the response data, p ji (t), and obtain the fitted probabilities, πg j (t), for each condition j. 2. Compute the errors for each trial i: ˆ ji (t) = p ji (t) − πg j (t), i = 1, . . ., I. 3. Estimate γ and σ 2 of the AR(1) model for the errors: j (t) = γ j (t − 1) + a t (as described in Section 4).Build the bootstrap samples under the null: 4. Fit the penalized general regression model under the null hypothesis to the response data, p ji (t), and obtain the fitted probabilities, π0 j (t), for each condition j.

Build bootstrap errors
) and z t ∼ N (0, σ2 ). 6. Compute the bootstrap data by adding the bootstrapped errors to the null model: p * ji (t) = π0 j (t) + * ji (t). 7. Fit the regression model from Step 1 to the bootstrap data to obtain the bootstrap synchrony curve π * j (t).8. Repeat Steps 5-7, M times to obtain M averaged bootstrap curves for each condition: π * 1 j (t), . . ., π * M j (t).For the first hypothesis test, the model fitted in Step 4 is the one that corresponds to the null hypothesis of constant synchrony before stimulation, the baseline model: logit(π pre (t)) = β 0 and it is fitted using only the data from the period previous to stimulation.Once the bootstrap curves are obtained, compute the appropriate percentiles to build either a poitwise or a uniform critical region around the baseline model for each condition.
For the second hypothesis test we need a further step.The null hypothesis H 2 0 states that π 1 (t) = π 2 (t) or, equivalently, that π 1 (t) − π 2 (t) = 0. So, in this case, we want to build critical bands for the difference of the synchrony curves.First of all, in Step 4, we need to fit the model which represents this hypothesis, this is Model (3) described in Section 4.Then, we will subtract the bootstrap curves to obtain: After this additional step we will compute the appropriate percentiles to build the critical region (pointwise or uniform) around zero (difference under the null hypothesis) and observe where the estimated difference curve πdif (t) = π1 (t)− π2 (t) falls inside that region, in case of searching pointwise testing.7.1.Pointwise and uniform critical regions.In order to assess whether at a given time point the observed CSM is different from the one observed at baseline at a level α, pointwise critical intervals can be built.To this aim, at each time point, we can choose the (1 − α/2) and α/2 percentiles of the bootstrap values: π * 1 j (t), . . ., π * M j (t).These marginal critical intervals have approximately a (1 − α) coverage probability at each time point separately.However, when considered as a uniform critical band, they do not provide an α-critical band for the null hypothesis: the curve is different from baseline over a whole period of time after stimulation because of the multiple range testing problem.This problem arises from the fact that, although at each time point the interval defined by the (1 − α/2) and α/2 percentiles include approximately the (1 − α)% of the bootstrap curves, when considered as a whole, the actual coverage will be considerably smaller.That is, the probability that the whole CSM curve is included in the band formed by the individual intervals is smaller than (1 − α).A posible approach would be to use a Bonferroni correction [2] which proposes to build the marginal intervals with a new significance level given by α/s, where s is the total number of tests to compute.Nevertheless, this approach is known to be very conservative.To get more accurate uniform critical bands we use an iterative procedure proposed by [4].The method consists on starting with the critical regions given by the (1 − α) marginal intervals and the conservative critical band given by the Bonferroni method and iterate until finding an approximate (1 − α)% coverage critical band up to an approximating error given by τ (for example τ = α/10).The algorithm to find this critical band, as presented in [4], is: high .These proportions satisfy: Otherwise increase k by one unit and repeat Steps 2-5.

Results.
When the firing rate is high, or a large amount of trials are available, both parameters v and δ can be specified small.On the other hand, if there are only a few trials or the firing rates are very low, wider windows need to be chosen to ensure that spikes are present in the intervals.How to choose the quantities v and δ is a matter of discussion and the effect of these selections will be exemplified in Section 9. We will see that the selection of v is not really relevant if working with GAMs.For our particular example, we use the smallest values that will allow computations of the CSM. Figure 3 shows an example carried out with the real data.In this figure, we observe the amount of zeros (no presence of CNSIs smaller than δ) encountered while computing the CSM on a 150 seconds trial of a real neuron.We have repeated computations for 45 pairs (v,δ): v = 3, 4, 4.5, 5, 5.5 seconds and δ = 0.005, 0.01, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05, 0.1 seconds.Each line corresponds to one value of v.We see that the amount of zeros for small windows together with small δ is large.In what remains, we will use v = 5 s and δ = 0.05 s.In our example, the main interest of studying synchrony between neurons lies in its time evolution.In this section we present an example in which synchrony is aimed to be observed to find differences between experimental conditions.Neurons are under spontaneous activity and therefore their firing rate is very low.As the stimuli are applied in a specific time point and the rest of the recordings correspond to spontaneous activity, our intention is to determine in what periods of time the CSM profiles differ.
The models described in Section 4 were fitted and the results are presented here.In this context, we have 2 experimental conditions which stand for stimulation in bs and bf.In Table 1 the AIC values for each model can be observed and, based on them, we can say that Model 1 is a better one.On the other hand, the intercept for stimulus 1, β1 , is not significant.
In Figure 4, we can observe that the curve that accounts for the difference between experimental conditions in the GAM, f 1 (t), is different from zero.This means that there are differences between the conditions and different smooth terms are needed for each stimulus.This hypothesis is tested statistically using the bootstrap procedure discussed in Section 7 and the results will be discussed shortly.Before that, in    Figure 7 shows the results for the second hypothesis test, using M = 500 bootstrap replicates, and although the critical band for zero (expected synchrony difference under the null hypothesis) is quite wide, we can still find a time period where the observed difference lies outside the band, meaning that the differences in synchrony observed between stimuli are significant.Finally, Figure 8 shows the uniform 95%-acceptance bands for a period of time after stimulation where the hypothesis of difference is relevant.The uniform acceptance band was constructed, using the procedure described in Section 7.1 for each of the hypotheses, as described in Section 5. A total of 10,000 bootstrap curves were resampled to construct each band.The final levels used to get the approximate 95% coverage of the bootstrap curves for the baseline tests were: 0.0074 for bs and 0.0073 for bf, achieving 95.04% and 95.16% final coverage, respectively.For the hypothesis test regarding the possible difference between the effect of the stimuli, the final level used was 0.0075, achieving a coverage of 94.91%.In the case of bs we can confirm that, overall, the synchrony curve after the stimulus onset is different from baseline synchrony, although we cannot assess this for the bf stimulus.We can also assess that the difference between synchrony curves after stimulation is significantly different from zero. .CSM profiles (solid lines) in an interval after stimulation with uniform 95%-acceptance bands (dashed lines) for the null hypotheses of no difference from baseline after bs stimulation (upper panel) and after bf stimulation (middle panel).Also, for the null hypothesis of no differential effect of the stimuli (bottom panel) Next, we show the results for the comparison of the observed synchrony to the expected one by chance.Figures 9 and 10 show these results for the real data.The first one shows the calculations from the analytical expression in 4 and the second one the results of the simulations.In both cases we can observe that the expected synchrony under independence is smaller than the observed one.There are short periods of time where the synchrony observed by chance gets close to the observed one, for example in both figures for the bs stimulus right after the stimulation has taken place.This uprise in the expected synchrony is due to the increase of the firing rates so the real synchrony might not have increased but actually have decreased.According to our working hypothesis we would expect the synchrony to decrease when the stimulus is applied.9. Simulation study.In this section we study the performance of CSM to capture neural synchrony on simulated data.We describe the effect of the tuning parameters, v and δ, on the measure.Also, we present a simulation study to evaluate the performance of the bootstrap test.
The data used for this aim was simulated using Poisson processes with controlled association.To generate a pair of simultaneous spike trains with a given firing rate λ, we first generated a master Poisson sequence of spike times with firing rate λ M .Then, the first of the two desired spike trains was constructed as follows.Given a real number q ∈ [0, 1], for each spike time of the master process, a random number, q, was drawn with uniform probability from [0, 1] and compared to q.If q < q, the spike time was inherited by the spike train.The spike times were then jittered by adding a small amount of noise drawn from a uniform distribution in ( −1 10λ , 1 10λ ).The whole process was repeated to construct the second spike train.Note that the probability of both trains inheriting the same spike is q 2 .Thus, the association of the two trains is controlled by the probability q.The final firing rate of the generated spike trains is also of interest.The firing rate λ is directly associated with the firing rate of the master process: λ = λ M q.Therefore, we chose λ M = λ q so that the simultaneous spike trains have approximately the desired firing rate.9.1.Performance of the method.Figure 11 shows the performance of CSM for different values of q.A constant firing rate of 4 Hz was used to simulate the spike trains and three choices of q were used.Spike activity was simulated for 300 seconds, 100 seconds for each value of q.The CSM was computed for two different choices of δ and v.The top left panel shows the CSM of one pair of simulated data.The effect of the change in the association parameter is very clear.In this case, we observe that δ and v have a very differential effect on the measure.Larger values of v give a smoother curve than smaller values of v, while the effect of δ is of another nature.As expected, larger values of δ values increase the value of CSM, as the proportion of CNSIs smaller than δ over the whole amount of CNSIs is larger.The top right panel shows the CSM averaged over 200 pairs of simulated spike trains.The effect of the change in the association parameter is even clearer in this figure.The effect of δ can still be observed while the effect of v is clearly diminished when averaging.The bottom panels show the fitted GAM for the CSM for one pair (bottom left) and averaged over 200 pairs (bottom right).It is very clear in this case that the choice of v does not influence the outcome.As a matter of fact, the curves corresponding to the same δ value but different v value (green-black, redblue) are completely confounded.Nevertheless, δ still does influence the measure considerably.Therefore, the choice of δ is very relevant and it needs to be done carefully, taking into account the firing rates of the spike trains and, most of all, the type of study that one wants to carry out.An important issue is to note that Hz and 50 Hz) firing rates and step-wise decreasing association parameter q (q = 0.99, 0.5 and 0.0001).
The procedure was repeated 1000 times and the proportion of rejections of the null was computed for each time point.Results are shown in Figure 13.It can be observed that, when the difference in the association parameter, q, is large, the test performance is really good.On the left hand side of each graph the proportion of rejections is almost zero, while in the right half, the proportion grows to one.In a time interval around 50 s the proportion continuously grows from zero to one as expected due to the use of sliding windows.When the association parameter used in both halves of each panel is small, the performance is not as good.The proportion of false positives increases as well as the proportion of false negatives diminishes.10.Discussion.A method to measure neural synchrony dynamics, specially under low firing rates scenarios, has been proposed in this work.The method is based on the computation of times between nearest spikes of two spike trains.This is a flexible method as it allows the researcher to decide how close two spikes have to be, to be considered synchronous.The method uses the information of time windows to estimate synchrony at a given time point.Although temporal resolution is lost with the use of this kernel-like method, it allows to study synchrony even with very few -or even only one-trials.Generalized additive models are proposed for the curves, giving flexibility to their shape by the use of nonparametric functions.Nevertheless, the dependency between consecutive time points is taken into account by modeling

Figure 1 .
Figure 1.Firing rates of two trials of two simultaneously recorded neurons.A trial in which the brainstem was stimulated (top panel) and a trial in which the basal forebrain was stimulated (bottom panel).Stimulation time is represented by the vertical dotted line (50 s).The spike times are depicted under each plot.The top row of spike times corresponds to the solid-line firing rate function and the bottom one to the dashed-line function

Figure 2 .
Figure 2. Joint peristimulus time histogram computed across three trials corresponding to the bs stimulus for different bin widths: 1 ms (top panel) and 100 ms (middle panel).CSM estimated with the same three trials (bottom panel)

3 . 4 .
Use the bootstrap resamples to build the marginal critical intervals with α Compute the proportion of simulated curves included in each of the critical bands: ρ

Figure 3 .
Figure 3. Number of zeros obtained in the calculation of CSM when using different δ (x-axis) and different v for the time windows ((t − v, t + v]): v = 3 s (black), v = 4 s (red), v = 4.5 s (green), v = 5 s (blue) and v = 5.5 s (cyan)

Figure 5
Figure5we observe the bootstrap confidence bands built for each of the estimators using M = 500 bootstrap replicates.

Figure 4 .Figure 5 .
Figure 4. Smooth terms for the GAM (solid lines) for CSM with confidence intervals (dashed lines).Mean curve for both experimental conditions f 0 (t) (upper panel) and smooth curve that accounts for the difference between experimental conditions, f 1 (t), (bottom panel)

Figure 6 .
Figure 6.CSM profile (thick solid lines) for two stimuli averaged over three trials: bs (upper panel) and bf (bottom panel).Baseline CSM estimated from the pre-stimuli time period (dotted lines) with bootstrap critical bands (thin solid lines)

Figure 7 .
Figure 7. Difference of CSM profile (thick solid line) for two stimuli averaged over three trials and bootstrap confidence bands (thin solid lines) for the null hypothesis of difference equal to zero (dotted line)

Figure 8
Figure 8. CSM profiles (solid lines) in an interval after stimulation with uniform 95%-acceptance bands (dashed lines) for the null hypotheses of no difference from baseline after bs stimulation (upper panel) and after bf stimulation (middle panel).Also, for the null hypothesis of no differential effect of the stimuli (bottom panel)

Figure 9 .
Figure 9. CSM profile (black solid lines) for bs stimulation (upper panel) and bf stimulation (bottom panel) compared with the theoretical expected synchrony (red solid lines).The moment of stimulation is indicated with vertical dotted lines.Also, the confidence bands for the CSM estimator are shown in this figure (dashed lines)

Figure 10 .
Figure 10.CSM profile (black solid lines) for bs stimulation (upper panel) and bf stimulation (bottom panel) compared with the simulated expected synchrony (red lines).The moment of stimulation is indicated with vertical dotted lines.Also, the confidence bands for the CSM estimator are shown in this figure (dashed lines)

Figure 12 .
Figure 12.Generalized additive model for CSM computed with different values of δ (0.05 s and 0.005 s) for pairs of spike trains simulated with different constant (4 Hz, 12 Hz and 50 Hz) firing rates and step-wise decreasing association parameter q (q = 0.99, 0.5 and 0.0001).

Table 1 .
AIC values for the two models